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INTRODUCTION 

The  concept  of  the  boxmdary  layer  was  first  introduced  by  Prandtl 
(Ref.  35)  so  that  approximate  solutions  could  be  obtained  to  the  Navier- 
Stokes  equations  which  describe  the  behavior  of  viscous  flow  around 
moving  bodies.  This  concept  allows  one  to  divide  the  flow  field  around 
a  body  into  two  parts.  In  the  external  part  the  effect  of  the  viscosity 
of  the  fluid  is  neglected  and  the  Navier-Stokes  equations  reduce  to  the 
Euler  equations.  In  the  other  part  of  the  flow  field,  which  is  near 
the  body,  the  viscosity  has  a  strong  influence;  but  other  terms  in  the 
Navier-Stokes  equations  can  be  neglected  to  give  the  classical  boundary- 
layer  equations.  As  a  result  of  the  simplification  of  the  Navier-Stokes 
equations  to  the  boundary- layer  equations,  the  pressure  change  normal 
to  the  wall  across  the  layer  is  negligibly  small.  Therefore  the  pres¬ 
sure  in  the  boundary- layer  equations  is  replaced  by  the  pressure  at  the 
wall  as  determined  from  the  external  inviscid  fluid  problem.  With  the 
tangency  boundary  condition  applied  at  the  surface  of  the  body,  the 
solution  of  the  Euler  equations  therefore  gives  the  pressure  distribution 

J 

on  the  body,  and  hence  the  pressure  gradient  is  known  in  the  boundary-)- 
layer  problem. 

The  presence  of  a  bo\indary  layer  on  a  body  effectively  increases 
the  thickness  of  the  body  by  an  amount  equal  to  the  displacement  thick¬ 
ness  as  far  as  the  external  flow  is  concerned.  Since  the  layer  thick¬ 
ness  is  small,  the  pressure  field  of  the  thickened  body  and  of  the  body 
proper  are  practically  equeLL  in  subsonic  and  low  supersonic  flow.  How¬ 
ever,  at  hypersonic  speeds  the  boundary  layer  cannot  be  neglected  in 
solving  the  Euler  equations  for  the  external  flow  about  slender  bodies. 
Since  the  displacement  thickness  is  a  function  of  the  press\ire  distribu¬ 
tion  along  the  body  and  the  pressure  distribution  is  a  function  of  the 
effective  body  shape,  there  is  an  interaction  between  the  boiindary- layer 
problem  and  the  external  flow  problem.  Therefore,  the  E\iler  equations 
and  the  bo\indary- layer  equations  must  be  solved  simultaneously. 

If  a  formal  approach  is  applied  to  approximate  the  Navier-Stokes 
equations  (technique  of  inner  and  outer  expansidns  as  applied  by 

-  1  - 


by  Van  Dyke  (Ref.  to  a  blunt  body)  the  classical  boundary-layer 
equations  are  the  first  approximation  to  the  Navier-Stokes  equations 
for  the  boundary- layer  flow.  The  second  approximation  to  the  Wavier- 
Stokes  equations  for  the  bo\mdary- layer  flow  introduces  the  effect  of 
the  displacement  thickness.  Also  the  following  second-order  effects  are 
introduced:  longitudinal  curvatiire^  transverse  curvature^  sljip^  tempera¬ 
ture  Jump  at  the  surface^  entropy  gradient^  and  total  enthalpy  gradient. 

In  this  paper  only  the  boundary- layer  flow  downstream  of  the  leading  edge 
of  a  sharp  flat  plate  will  be  considered.  For  a  flat  plate  the  longi¬ 
tudinal  and  transverse  curvature  effects  are  zero.  Since  slip  and  temp¬ 
erature  Jump  at  the  surface  are  mainly  important  near  the  leading  edge^ 
the  boundary- layer  flow  is  considered  only  after  a  small  distance  from 
the  leading  edge.  For  a  sharp  flat  plate  there  will  be  a  nearly  straight 
shock  wave  from  the  leading  edge  and  hence  the  entropy  gradient  will  be 
small.  The  stagnation  enthalpy  is  constant  across  the  shock  wave  and 
hence  the  total  enthalpy  gradient  is  zero.  Therefore,  the  predominant 
second-order  effect  will  be  the  displacement  thickness  and  this  will  be 
the  only  effect  considered  in  the  interaction  between  the  boundary- layer 
and  the  external  flow. 

Before  the  interaction  problem  can  be  solved  it  is  necessary  to  have 

« 

a  satisfactory  method  to  solve  the  boiondary-layer  equations  with  an  arbi¬ 
trary  pressure  distribution.  Since  similarity  methods  are  only  applicable 
for  rather  special  pressure  distributions  and  integral  methods  only  give 
approximate  results^  the  numerical  scheme  of  finite-differences  is  em¬ 
ployed.  Several  finite- difference  methods  for  the  compressible  laminar 
boundary  layer  existed  when  this  work  started.  Baxter  and  Flugge-Lotz 
(Ref.  l)  modified  the  boundary-layer  equations  using  the  Crocco  trans¬ 
formation  and  then  used  ein  explicit  finite-difference  scheme  to  solve 
the  resulting  eqmtlons.  This  method  is  not  completely  satisfactory  due 
to  the  small  step -sizes  required  for  ensuring  stability  and  convergence 
of  the  finite-difference  scheme.  Kramer  and  Lleberstein  (Ref.  21)  have 
solved  essentially  the  same  Crocco  transformed  equations  except  an  im¬ 
plicit  finite-difference  scheme  has  been  used  to  eliminate  stability 
problems.  Still  the  step -size  must  be  sufficiently  small  to  ensure 
convergence  of  the  n\americal  solution  to  the  exact  solution.  In  their 
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paper  there  are  no  comparisons  vlth  known  results  or  any  indications  of 
the  validity  of  the  method.  A  definite  disadvantage  of  methods  based 
on  Crocco's  equations  is  the  fact  that  the  methods  are  not  applicable 
when  the  boundary- layer  profiles  have  "overshoot."  Velocity  profiles 
with  overshoot  occur  for  certain  cases  of  heated  walls  with  favorable 
pressure  gradient  and  this  effect  becomes  even  more  important  for  boundary- 
layer  flows  with  helliam  injection.  Due  to  the  above  considerations^ 
Flugge-Lotz  and  Yu  (Ref.  ll)  investigated  an  explicit  finite-difference 
scheme  to  solve  the  boundary- layer  equations  for  the  physical  plane. 
However,  this  method  did  not  prove  completely  satisfactory,  especially 
at  high  Mach  numbers  and  with  a  heated  wall.  Under  these  conditions 
the  step-size  requirements  were  so  severe  that  it  was  impossible  to  ob¬ 
tain  stable  solutions. 

Since  the  available  numerical  solutions  of  the  boundary- layer  equa¬ 
tions  were  not  generally  adequate,  a  major  part  of  this  paper  is  con¬ 
cerned  with  developing  a  new  method.  A  satisfactory  method  will  be  one 
which  is  not  only  stable,  but  one  which  has  a  grid  size  such  that  compu¬ 
tation  time  is  reasonable.  Two  main  implicit  finite-difference  schemes 
are  developed.  Since  the  boundary- layer  equations  are  analogous  to  the 
heat  equation  as  far  as  stability  is  concerned,  an  implicit  scheme  using 
the  boundary- layer  equations  in  the  physical  plane  was  used  initially 
in  order  to  have  stability  as  the  heat  equation  indicates.  One  of  the 
better  ways  to  solve  the  heat  equation  is  the  Crank-Nicolson  method 
(Ref,  8)  which  has  a  smaller  truncation  error  than  the  implicit  method. 

This  method  actually  is  at  the  dividing  point  between  the  explicit  and 
implicit  schemes,  but  is  always  stable  for  the  heat  equation.  Due  to 
the  complexity  of  the  boundary- layer  equations,  it  is  impossible  to  as¬ 
certain  with  assurance  whether  a  Crank-Nicolson  type  of  scheme  would 
always  be  stable.  In  order  to  improve  the  grid  size  requirements,  and 
thus  have  a  practical  method,  the  second  difference  scheme  investigated 
is  of  the  Crank-Nicolson  type. 

Especially  for  hypersonic  flows  along  Insulated  walls,  the  velocity 
profiles  vary  almost  linearly  over  a  large  portion  of  the  layer  except 
at  the  outer  edge  where  there  is  a  large  change  in  the  velocity  gradient. 
For  accurate  results  a  particularly  small  grid  size  is  required  near  the 
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outer  edge.  The  immediate  conclusion  is  that  the  grid  size  should  vary 
with  distance  from  the  wall,  but  this  is  rather  inconvenient.  Therefore, 
the  coordinate  normal  to  the  wall  should  be  stretched  to  obtain  smoother 
profiles  across  the  boundary  layer  by  using  the  Howarth-Dorodnitsyn 
transformation.  This  was  done  in  the  second  scheme  investigated. 

In  both  of  the  implicit  finite-difference  schemes  developed  in  this 
report  the  derivatives  in  the  boundary- layer  equations  are  replaced  by 
difference  quotients  such  that  linear  difference  equations  are  obtained. 

If  the  usual  procedures  are  followed  for  replacing  the  partial  differen¬ 
tial  equations,  nonlinear  difference  equations  are  obtained.  Hence,  one 
has  the  extremely  involved  problem  of  solving  forty  to  sixty  simulteuieous 
nonlinear  algebraic  equations  at  each  step  along  the  wall.  In  order  to 
make  the  implicit  method  feasible,  the  difference  equations  have  been 
lineaLTized,  which  then  requires  the  solution  of  a  large  number  of  simul¬ 
taneous  linear  algebraic  equations.  These  equations  are  of  a  tridiagonal 
type  and  are  well  suited  for  solution  on  a  digital  computer.  A  wide  var¬ 
iety  of  boundary- layer  problems  has  been  solved  to  indicate  the  useful¬ 
ness  and  accuracy  of  the  two  methods  investigated.  The  problems  used  for 
the  verification  of  these  methods  either  have  exact  solutions  or  have  been 
solved  numerically  by  other  procedxxres. 

Since  the  above  finite-difference  schemes  have  been  developed, 
several  numerical  methods  to  solve  the  boundary- layer  equations  have 
appeared.  Wu  (Ref,  42)  has  used  an  explicit  finite- difference  scheme  to 
solve  the  boundary- layer  equations  in  the  physical  plauie.  The  Howarth- 
Dorodnitsyn  transformation  has  also  been  used  by  him  to  improve  the  sta¬ 
bility  requirements  for  the  compressible  boundary  layer.  However,  the 
transformed  equations  as  given  cannot  be  used  for  a  flow  with  a  pressure 
gradient;  this  fact  is  not  clearly  Indicated  in  the  report.  Another 
method  which  combines  the  Dorodnitsyn  integration  scheme  with  the 
Pohlhaiisen  approach,  has  been  developed  by  Pallone  (Ref.  34).  The  boun¬ 
dary  layer  is  divided  into  a  number  of  strips  (4  or  6)  parallel  to  the 
wall  and  then  the  boundary- layer  equations  are  integrated  from  the  wall 
to  the  various  strips  where  the  velocity  and  enthalpy  profiles  have  been 
approximated  by  a  polynomial.  This  method  reduces  the  partial  differen¬ 
tial  equations  of  the  boundary  layer  to  a  set  of  ordinary  first-order 
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differential  equations  with  the  coordinate  parallel  to  the  wall  as  the 
independent  variable.  This  method  looks  very  attractive  as  the  proce¬ 
dures  to  solve  ordinary  differential  equations  are  highly  developed. 
However^  one  would  question  whether  this  method  converges  to  the  exact 
solution  as  the  number  of  strips  is  increased,  since  a  polynomial  is  being 
fitted  through  a  large  number  of  points  to  represent  the  velocity  and  en¬ 
thalpy  profiles.  For  this  situation,  the  polynomial  can  be  greatly  dif¬ 
ferent  from  the  exact  profiles  between  the  fitting  points.  An  interest¬ 
ing  numerical  method  for  solving  the  incompressible  boundary- layer  flow 
has  been  presented  by  Manohar  (Ref.  29)  and  the  ideas  can  be  extended 
to  compressible  flow.  The  author  develops  from  the  boundary- layer 
momentum  equation  and  the  continuity  equation  a  third-order  nonlinear 
partial  differential  equation  and  then  replaces  the  derivatives  in  the 
direction  along  the  wall  by  difference  quotients.  This  results  in  a 
nonlinear  ordinary  differential  equation  which  must  be  solved  across 
the  boundary  layer  at  each  step  downstream.  An  iteration  process  is  re¬ 
quired  for  the  solution  since  two  boundary  conditions  are  given  at  the 
wall  and  one  at  the  outer  edge.  The  problem  of  iteration  would  become 
even  more  difficult  for  the  compressible  boundary- layer  equations  as 
there  would  be  two  boundary  conditions  at  the  outer  edge. 

A  large  nximber  of  people  have  investigated  the  displacement  thick¬ 
ness  or  pressiire  interaction  between  the  boundary-layer  and  the  external 
flow.  A  complete  discussion  of  this  problem  with  a  review  of  previous 
contributors  is  given  in  Hayes  and  Probstein  (Ref.  l4).  All  of  the  prev¬ 
ious  methods  of  solution  have  either  used  perturbation  or  approximate 
methods  to  solve  the  Prandtl  boundary- layer  equations  with  the  pressure 
gradient  in  these  equations  determined  from  the  effective  body  shape. 

In  this  paper  the  same  problem  is  investigated  using  the  new  numerical 
scheme  developed  for  solving  the  boundary-layer  equations.  The  tangent- 
wedge  formula  is  used  as  the  solution  to  the  external  flow  field  and  hence 
the  pressxrre  distribution  is  known  once  ^he  effective  shape  of  the  body 
(geometric  shape  plus  displacement  thickness)  is  given.  Starting  with 
the  Initial  profiles  across  the  boundary  layer,  the  boundary- layer  equa¬ 
tions  are  solved  with  the  pressure  at  the  next  grid  line  iterated  until 
the  assumed  pressure  equals  the  press\ire  determined  from  the  tangent- 
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wedge  form\J.a.  This  process  is  repeated  as  one  steps  downstream  until 
the  desired  distance  has  been  covered.  A  comparison  of  the  numerical 
results  is  made  with  the  "strong"  and  "weak"  interaction  theories  and 
experimental  results  when  such  are  available. 


CHAPTER  II 


MATHEMATICAL  DESCRIPTION  OF  THE  COMPRESSIBLE 
LAMINAR  BOUNDARY  LAYER 


In  this  chapter  the  equations  for  the  laminar  compressible  boxmdary 
layer  are  presented  along  with  the  necessary  boundary  conditions.  The 
equations  are  nondimens ionalized,  which  results  in  using  quantities  of 
the  same  order  of  magnitude  better  suited  for  computations .  The 
boundary- layer  equations  are  also  modified  using  the  Howarth-Dorodnitsyn 
transformation,  which  results  in  equations  advantageous  for  numerical 
computation  when  the  flow  is  hypersonic.  Since  the  boundary  conditions 
require  that  the  enthalpy  sind  velocity  be  known  at  the  outer  edge  of 
the  boundary  layer,  the  necessary  formulas  are  presented  to  calculate 
the  exterior  flow,  provided  the  pressure  distribution  is  known.  Finally, 
the  characteristic  values  of  the  boundary  layer  (shearing  stress,  heat 
transfer,  and  displacement  thickness)  are  defined. 

A.  The  System  of  Partial  Differential  Equations 

1.  Differential  Equations.  The  flow  of  a  compressible,  viscous, 
heat-conducting  fluid  is  mathematically  described  by  the  continuity, 
Navier-Stokes,  and  energy  equations  plus  the  equation  of  state,  a  heat- 
conductivity  law,  and  a  viscosity  law.  For  flows  at  large  Reynolds 
numbers,  Prandtl  (Ref.  35)  has  shown  how  the  Navier-Stokes  and  energy 
equations  can  be  simplified  to  the  boundary- layer  equations.  If  it  is 
assumed  the  Prandtl  number  and  specific  heat  of  the  fluid  are  constant, 
the  classical  boundary- layer  equations  for  steady,  two-dimensional, 
compressible  flow  are 


Continuity 

(p*u*)^^  +  (P*v*)y^^  =  0 

(2.1) 

<iPg* 

Momentum 

+  p*y*u*^  -  -  dx* 

y* 

(2.2) 
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Energy 


2  1 

p*u*i*^^  +  +  pj:  ^  (2-3) 

y 


Equation  of  State 


p*  =  p*  R*  T*  =  ^  p*i* 


(2.4) 


Viscosity  Law 


H*  =  f*(i*) 


(2.5) 


The  starred  quantities  have  dimensions  and  the  subscripts  indicate  par¬ 
tial  differentiation.  The  independent  variables  are  the  space  coordi¬ 
nates  X*  and  y*  which  are  parallel  and  perpendicular  to  the  wall, 
respectively.  The  dependent  variables  are  the  density  p*,  the  velocity 
components  u*  and  v*  which  are  parallel  and  perpendicxilar  to  the  wall, 
respectively,  enthalpy  i*  and  viscosity  p*.  Hence,  we  have  a  system  of 
five  simultaneous  equations  and  five  unknowns. 

The  two  original  momentxmi  eqmtlons  became,  after  the  boundary- 
layer  simplification,  equation  2.2  and  the  following  result 


=  0 


(2.6) 


Thus  the  pressure  in  a  direction  normal  to  the  boundary  layer  is  practi¬ 
cally  constant.  This  result  has  been  used  in  equations  (2.^)  and  (2.3) 
where  the  pressure  at  any  x*  has  been  set  equal  to  the  external  pressure 
p*.  Also,  from  the  eqmtlon  of  state  the  resiilt  is  obtained  that  p*!* 
is  a  constant,  or  the  following  relation  can  be  written 


(p*  ip/i^ 


(2.7) 


The  boundary- layer  equations  given  above  (2.1  -  2.5)  are  valid  for 
a  curved  wall  provided  there  are  no  large  variations  in  curvature,  and 
the  bovuidary- layer  thickness  is  smeill  compared  with  the  radius  of  curva¬ 
ture  of  the  wall. 

2.  Boundary  Conditions.  In  order  to  obtain  a  unique  solution  to 
the  partial  differential  equations  of  the  boundary  layer  it  is  necessary 
to  satisfy  the  boundary  conditions  of  the  problem  under  consideration. 
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These  conditions  are  shown  schematically  helow. 


The  velocity  and  enthalpy  at  the  edge  of  the  boundary  layer  are  deter¬ 
mined  from  the  shape  of  the  body  by  using  inviscid  flow  theoi^  and 
are  discussed  in  Section  E  of  this  chapter. 

The  conditions  of  no  slip  and  no  suction  or  blowing  are  taken  as 
the  boundary  conditions  at  the  wall.  Hence,  the  two  velocity  compo¬ 
nents  u*  and  V*  are  zero  at  the  wall.  The  third  boundary  condition  at 
the  wall  depends  upon  the  thermal  state  of  the  wall.  The  temperature 
of  the  wall  can  be  specified,  which  means  the  enthsLLpy  along  the  wall 
is  given.  An  alternate  condition  at  the  wall  is  to  specify  the  heat 
transfer  which  determines  .  Finally,  the  parabolic  charac¬ 

ter  of  the  differential  equations  requires  that  initial  velocity  and 
enthalpy  profiles  be  known  at  x*.  These  inltieLl  profiles  are  obtained 
from  analytical  results  which  are  applicable  to  the  problems  under 
consideration. 


In  the  classical  boundary- layer  theory  the  shape  of  the  body  is 
taken  as  the  geometriceil  body.  When  interaction  between  the  boundary- 
layer  and  inviscid  flow  is  considered,  the  shape  of  the  body  is  the 
effective  body  (see  Chapter  5) • 
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3.  Viscosity  Lav.  A  formiila  for  the  viscosity  of  a  gas  from  the 
standpoint  of  kinetic  theory  of  gases^is  Sutherland's  law,  which  is 


T*  +  S* 
o _ 

T*  +  S* 


(tVt*) 


3/2 


(2.8) 


The  constant  S*  is  tedcen  as  216  or  198.6  in  this  work  in  order  to  com¬ 
pare  results  obtained  with  other  reports  using  these  values.  Based 
on  newer  experiments  the  value  of  S*  is  198.6. 

Another  viscosity  law  which  is  used  in  many  analytical  solutions 
to  the  boundary- layer  equations  is  the  linear  law,  which  is  written  as 

nVti*  =  ctVt*  (2.9) 


The  constant  C  is  usually  determined  by  matching  the  viscosity  at  the 
wall  as  given  by  the  linear  law  and  Sutherland's  law.  This  results 
in  the  following  value  for  C: 


C 


T*  +  S* 
e 

e  T*  +  S* 
w 


Using  the  relation 


i*  c*  T* 
P 


(2.10) 


(2.11) 


gives  the  above  viscosity  laws  as  a  function  of  enthalpy. 


B.  Nondimens lonal  Form  of  the  Differential  Equations  and  Boundary 
Conditions 


1.  Nondlmenslonal  Queintities.  It  is  advantageous  to  write  the 
differential  equations  and  other  relations  of  interest  in  nondlmenslonal 
form  before  numerical  computations  are  performed.  The  following  non- 
dimensloned.  quantities  are  defined  so  that  all  quantities  are  of  the 
same  magnitude. 


Sutherland's  Hypothesis  assumes  that  the  moleciiles  are  hard  spheres 
with  a  weak  attraction  between  them  and  this  force  falls  off  rather  rapid¬ 
ly  with  distance.  The  resulting  viscosity  law  from  kinetic  theory  has 
the  constant  S*  which  is  a  suitable  empirical  value.  This  law  is  valid 
over  a  certain  temperature  range  (Bef .  19) • 
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a.  Variables  appearing  directly  in  the  differential  equations; 
u  =  u*/c* 


V  =  V* 


1  = 


P  =  P*/p* 


(2.12) 


with 


•p  =  pV(p*  P*  ) 


X  =  x*/l* 


y  =  y*7Re*/L* 


i^  =  l/(7-l) 


^0  =  1 


b.  Quantities  describing  the  behavior  of  the  boundary- layer 


flow: 


T  =  T* 


'=f'> 


q  =  q*/^Re*/(p*  cf) 


8=6*  S*/c*'^ 
p  '  o 

■^o 

r 

T  =  c  * 

p  '  o 


(2.13) 


h  =  h*Aj^lip*  6*  Cp*) 


(2.13  cont'd.  on  next  page) 


11 


(2.13  cont^d) 


— f 


k  =  kV(c  *  a*) 

c  =  c*/c  *  =  1 
P  p'  Pq 

Pr  =  Cp  ^i/k  =  c*  ^*/k*  =  Pr* 


R  =  R*/c  * 

P 

In  the  above  quantities  the  subscript  "o"  refers  to  the  conditions 
at  the  stagnation  point  or  reservoir  of  the  exterior  inviscid  flow.  The 
length  L*  is  some  characteristic  distance  of  the  problem  and  the  Reynolds 
number.  Re*,  is  defined  by 

Re*  =  c*  L*  p*/u*  (2.14) 

2.  Bo \jndary- Layer  Equations.  Since  the  viscosity  is  a  function 
of  enthalpy,  the  viscosity  in  nondimens ional  form  can  be  written  as 

^  =  f(i)  (2.15) 


The  derivative  of  the  viscosity  can  be  written  as 


y 


i 

y 


(2.16) 


Using  the  nondimens ional  quantities  as  given  in  (2.12),  we  can  write 
the  boimdary- layer  equations  (2.1  -  2.3)  as 

Continuity  (pu)  +  (pv)  =  0  (2.17) 

A  y 

Momentum  puu  +  pvu  =  -  p'  +  fu  +  f  u  i  (2.l8) 

x^y  ^e  yy  yy  '''' 

Energy  pui  +  pvi  =  up'  +  fu^  +  .^l  (2.19) 

X  ^  y  *^e  y  Pr  yy  Pr  y  \ 

Equation  of  State  p  =  p^i^/i  =  p/i  (2.20) 
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Viscosity  Law 


— t 


(Linear) 


t  =  C^± 


f  =  C, 


where  C.  =  Cu,  /i 
1  e'  e 


1  1  +  s/t 

w  '  '■ 


i„  1  +  s/t 

°  l“  i  /i  +  S/T 
e  w'  e  '  e 


(Sutherland) 


f 


V. 


1  +  s/t 


Vlo  +  S/T„ 


°  (1/i  )3/2 


'l/l„  3  S/T„ 

2i(i/i^+  S/T^) 


(2.22) 


3.  Boundary  Conditions.  The  boundary  conditions  become  the  fol¬ 
lowing  when  the  nondimens ional  quantities  are  used: 


a.  Outer  Edge  of  Boundary  Layer 


y  =  ye  ,  x>Xi  : 

u(x,yg)  =  Ug(x) 

i(x,yg)  =  lg(x) 

(2.23) 

Wall 

o 

II 

o 

X 

y  =  0  ,  X  >  x^  : 

v(x,o)  =  0 
i(x,o)  =  i^ 

iy(x,o)  .  (ly)^ 

(2.24) 

Initial  Profiles 

ii(Xi,y)  =  u^(y) 

:  v(x^,y)  =  v^(y) 

i(x^,y)  =  ij^(y) 

(2.25) 
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C.  Howarth  -  Dorodnltsyn  Transformed  Boundary-Layer  Equations 


For  some  numerical  calculations  it  is  advantageous  to  stretch  the 
coordinate  normal  to  the  wall.  The  appropriate  stretching  is  accomplished 
by  using  the  Howarth  -  Dorodnitsyn  transformation 


I  =  X 


n  = 


Pdy 


(2.26) 


In  order  to  transform  the  boundary-layer  equations  (2,17  -  2.19)  "ttie 
following  relations  between  the  differential  operators  in  the  old  and 
new  coordinates  are  usedc 


^  A  ^ 

Also,  a  new  velocity  is  introduced  and  is  defined  as 

V  =  u  +  pv  (2.28) 

The  botindary- layer  equations  (2.I7  -  2.19)  become  the  following  in 
the  transformed  pleme: 


=  0 

(2.29) 

UU|  +  V  u^  = 

=  -  p'/p  +  F  u^„  +  F'  i 
e'  Vn  T1 

U 

TJ 

(2.30) 

uig.Vi^. 

0  p 

up '  /p  +  F  u^  +  —  i  + 
“^e/^  ^  Pr  TjTi 

F'  2 

Pr  ^T) 

(2.31) 

F 

=  pu  =  pgig  (p/i) 

F' 

'  §  ■  31 

(2.32) 
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Of  course,  the  equation  of  state  remains  the  same,  but  the  viscosity- 
law  gives  the  following  relations; 


(Linear  Law) 


(Sutherland  Law) 


r. 


F  = 


p  i  C, 
^e  e  1 


=  o 


(2-33) 


(  1  +  S/T^  \ 

■  ip  vi„  *  s/T„y 

/S/T,  -  1/lA 

■  l^S/T^  +  1/1 J 


1/2 


(2.34) 


In  order  to  determine  the  physical  distance  normal  to  the  wall,  y,  the 
trsinsformation  (2.26)  is  used  to  obtain  the  following; 


y  = 


f 


(l/p)dTl  = 


(1/Pe) 


f 


(i/le) 


dT] 


(2.35) 


D.  Exterior  Flow 

The  following  quantities  will  be  given  to  describe  the  exterior 
flow  and  fluid  properties; 


M  ,  T  ,  7,  Pr,  and  S. 

oo'  oo'  ' 


Since  we  know  the  above  quantities,  the  following  formulas  for  an 
isentropic,  perfect  gas  can  be  evaluated  independent  of  the  shape  of 
the  wall  on  which  the  boundary  layer  is  located. 


u  =  M  1  +^~iiF 
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The  next  exterior  flow  quantity  that  must  he  determined  is  the  pressure 
and  pressure  gradient  along  the  boundary  layer.  The  method  to  be  used 
to  determine  this  pressure  depends  upon  the  problem  under  consideration; 
hence^  it  will  be  discussed  along  wlthtiie  problems  considered  in  later 
sections . 

Now  the  following  formula  for  an  isentropic,  perfect  gas  can  be 
evaluated  using  the  previously  obtained  results  above. 
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It  should  be  noticed  that  the  contribution  of  v^  in  equation  (2.4l) 
has  been  neglected, which  is  consistent  with  the  boundary- layer  equations. 


E.  Boundary- Layer  Parameters 

1.  Shearing  Stress.  Since  the  viscous  drag  of  a  body  moving  throiigh 
a  fluid  is  dependent  upon  the  shearing  stress  at  the  wall,  this  is  one 
of  the  boiindary- layer  characteristics  of  Interest.  The  shearing  stress 
at  the  wall  in  the  physiceil  plane  is 
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and  in  the  transformed  plane  is 


\  >v 

w 

Rather  than  the  above  parameter,  usmlly  the  skin  friction  and  Reynolds 
number  8J?e  combined  to  give  one  of  the  following  characteristics  to 
describe  the  boundary- layer  shear  stress  at  the  wall; 


where 


(2.47) 
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(local  skin -friction  coefficient) 


(local  Reynolds  number) 


In  the  case  of  interaction  between  the  boundeury- layer  and  the  external 
flow,  the  following  parameter  is  used  to  describe  the  skin  friction 
and  is  used  for  ansilytical  resialts  in  Ref.  l4. 
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2.  Heat  Transfer.  In  aerodynamic  heating  analysis  it  is  necessary 
to  know  the  heat  transfer  between  the  boundary  layer  and  weLLl.  The  heat 
transfer  at  the  wall  is  related  to  the  enthalpy  by  the  following  formula: 


S.  =  -  I"  ^^-^9) 

p  w  w 

Again,  there  are  more  conventional  parameters  to  describe  the  heat 
transfer  at  the  wall  and  one  is 


Pe  ^^e  % 
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where 


(the  Stanton  number) 

(local  heat  transfer 
coefficient) 


i  =  i  (1  +  ^  r  W^) 

The  adiabatic  enthalpy,  i^^,  is  the  value  of  the  enthalpy  at  the  wall 
when  there  is  no  heat  transfer  (insulated  wall).  In  order  to  obtain 
the  adiabatic  enthalpy,  it  is  necessary  to  know  the  recovery  factor  r. 
For  the  laminar  flow  along  a  flat  plate  at  constant  temperature  the  re¬ 
covery  factor  is  taken  as  0.845  when  Pr  =  0.72.  For  some  other  types 
of  flow  the  recovery  factor  is  unknown.  In  these  cases  the  Stanton  num¬ 
ber  can  be  based  upon  the  temperature  at  the  edge  of  the  boundary  layer, 
Tg,  and  then  all  the  adiabatic  enthalpies,  i^^,  are  replaced  by  i^  in 
the  above  relations . 

For  the  problem  of  interaction  between  the  boundary- layer  and  the 
external  flow,  the  following  parameter  is  used  to  describe  the  heat 
transfer  and  is  used  for  experimental  results  in  Ref.  (13 )• 
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3.  Displacement  Thickness.  For  problems  in  which  interaction  be¬ 
tween  the  boundary- layer  and  inviscid  flow  is  important,  it  is  neces¬ 
sary  to  have  the  displacement  thickness  to  determine  the  effective  body. 
The  boundary-layer  displacement  thickness  is  defined  as 


(2.52) 
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A  nondimens lonal  displacement  thickness  is  used  in  calculations  as  de¬ 
fined  below; 
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In  the  transformed  plane  the  nondimens ional  displacement  thickness  be¬ 
comes 
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CHAPTER  III 


NUMERICAL  SOLUTION  OF  THE  BOUHPARY- LAYER 
EQUATIONS  WITHOUT  INTERACTION 

Since  the  partial  differential  equations  governing  the  flow  in  the 
boundary  layer  eore  of  parabolic  type,  they  can  be  solved  stepwise  down¬ 
stream  starting  with  initial  velocity  and  enthalpy  profiles  and  using  the 
boundary  conditions.  When  the  finite-difference  scheme  is  employed,,  the 
derivatives  in  the  partial  differential  equations  are  replaced  by  dif¬ 
ference  quotients,  and  this  results  in  difference  equations.  To  con¬ 
struct  the  difference  quotients  and  equations  it  is  necessary  to  divide 
the  flow  field  into  a  grid  or  mesh.  First,  the  velocities  and  enthaply 
across  the  boundary  layer  at  the  grid  points  which  are  a  small  distemce 
downstream  from  the  initial  profiles  are  computed  using  the  velocities 
and  enthalpy  given  by  the  initial  profiles.  When  the  difference  equa¬ 
tions  are  solved  in  succeeding  small  steps  downstream,  the  flow  in  a 
region  of  the  boundary  layer  can  be  determined. 

There  are  nximerous  ways  of  constructing  difference  quotients,  but 
there  are  two  general  classifications:  explicit  and  implicit.  When 
the  explicit  scheme  is  used,  the  resulting  finite- difference  equations 
can  be  solved  for  the  unknown  quantities  at  one  grid  point  at  a  time. 
However,  to  obtain  stable  numerical  solutions,  the  step-size  downstream 
must  be  less  than  some  critical  value.  With  the  implicit  scheme  the  dif- 
:ference  equations  for  the  ixnknown  quantities  at  a  row  of  grid  points  are 
solved  simultaneously.  The  downstream  step-size  for  the  implicit  method 
requires  no  restrictions  to  ensure  stability;  however,  the  step-size 
must  be  sufficiently  smeill  for  the  numerical  solution  to  closely  approxi¬ 
mate  the  exact  solution.  Obviously,  the  implicit  scheme  is  mathematically 
more  complicated;  however,  whentiie  difference  equations  are  linear,  a  di¬ 
rect  method  of  solution  (an  algorithm)  can  be  evolved.  Due  to  the  repe- 
titioxis  form  of  the  explicit  and  implicit  finite- difference  calculations, 
the  solution  is  well  suited  for  digital  computers.  The  best  method  for 
a  computer  is  determined  from  the  consideration  of  computation  time 
needed  rather  than  the  mathematical  complication  involved  or  the  step- size 
required.  The  implicit  methods  usually  require  more  computation  time  per 
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step,  tut  allow  larger  step -sizes  downstream  than  the  explicit  methods; 
therefore,  total  computation  time  by  either  method  can  be  a  minimum  for 
a  particular  problem. 

The  essential  points  of  the  explicit  method  of  Wu  (Ref.  42)  are 
presented  in  this  chapter  for  convenience  of  comparison  with  the  impli¬ 
cit  methods.  Also  the  equations  are  presented  in  a  slightly  different 
form  which  is  valid  for  boundary- layer  flow  with  a  pressvire  gradient. 

The  equations  given  by  Wu  are  only  valid  for  zero  press\ire  gradient  since 
the  density  at  the  edge  of  the  boundary  layer  has  been  assumed  a  con¬ 
stant  while  applying  the  Howarth-Dorodnitsyn  transformation.  The  ex¬ 
plicit  method  is  presented  only  for  the  solution  of  the  boimdary- layer 
equations  in  the  transformed  plane  (after  the  Howarth-Dorodnitsyn  trans¬ 
formation  has  been  used).  In  the  explicit  scheme  the  derivatives  normal 
to  the  wall  are  replaced  using  quantities  from  the  known  profiles.  There¬ 
fore,  the  momentum  equation  gives  the  velocity,  u,  and  the  energy  equa¬ 
tion  gives  the  enthalpy,  i,  directly  at  each  grid  point  (except  at  the 
edges)  one  step  downstream  from  the  known  profiles.  The  quantities  at 
the  edges  of  the  boundary  layer  are  known  from  the  boundary  conditions. 

The  velocity,  V,  is  determined  from  the  continuity  equation  and  then 
the  above  procedure  is  repeated  at  succeeding  steps  downstream.  Wu 
has  Indicated  that  the  continuity  equation  in  his  method  is  treated  in 
a  particular  manner.  The  derivative  u^  is  replaced  with  a  backward 
difference  quotient  and  the  resulting  equation  is  integrated  using  the 
trapezoidal  rule. 

Two  implicit  finite-difference  schemes  are  presented  in  this 
chapter  for  the  solution  of  the  nonlinear  boundetry- layer  equations. 

The  first  scheme  is  the  usual  implicit  one  which  uses  backward  difference 
quotients  for  the  derivatives  parsLllel  to  the  wall  and  is  called  Method  I. 
The  second  scheme  is  of  the  Crank-Nicolson  type  which  usee  difference 
quotients  of  a  type  to  reduce  the  truncation  error  and  is  called  Method  II. 
Both  Methods  I  and  II  replace  the  derivatives  normal  to  the  wall  in  such 
a  manner  that  unknown  quantities  eu:*e  introduced.  Also,  the  products  of 
the  derivatives  are  replaced  with  linear  difference  quotients  to  give 
the  possibility  of  linear  difference  equations.  However,  to  obtain 
linear  difference  equations,  it  is  necessary  to  linearize  certain  terms 
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obtained  from  the  partial  differential  equations.  In  the  process  of 
linearizing  the  momentum  and  energy  equations,  only  the  velocity,  .u, 
and  enthalpy,  i,  appear  as  unknowns  in  the  resulting  difference  equa¬ 
tions.  In  both  methods  these  difference  equations  are  of  the  same  form 
with  only  the  coefficients  depending  on  which  method  is  being  used. 
Therefore,  the  manner  of  solution  is  identical  in  both  cases  with  the 
appropriate  coefficients  used  for  the  two  methods.  Since  at  the  first 
grid  point  away  from  the  wall  the  two  difference  equations  have  six 
unknowns  and  the  boundary  conditions  at  the  wall  can  eliminate  two 
unknowns,  there  are  more  unknowns  than  equations.  If  the  two  differ¬ 
ence  equations  are  introduced  for  the  next  grid  point,  there  will  still 
be  more  unknowns  than  equations.  There  will  always  be  more  unknowns 
than  equations  as  all  the  difference  equations  are  added  for  the  grid 
points  across  the  boxindary  layer  until  the  outer  edge  is  reached. 

Using  the  boundary  conditions  at  the  outer  edge  with  those  at  the  wall, 
one  obtains  a  set  of  difference  equations  with  the  same  number  of  un¬ 
knowns  as  equations.  This  set  of  difference  equations  is  actually  a 
rather  special  system  of  simultaneous  linesu:  algebraic  equations.  Due 
to  this  special  form  an  efficient  method  of  solution  for  computers  is 
available;  and  using  this  method,  the  difference  equations  are  replaced 
by  suitable  relations  useful  for  computation  in  this  chapter.  This 
method  essentially  introduces  six  new  quantities  which  can  be  determined 
directly  by  using  the  boundary  conditions  at  the  wall  and  then  stepping 
across  the  layer.  Then,  using  the  boundary  conditions  at  the  outer  edge 
and  these  six  new  quantities,  one  determines  the  velocity,  u,  and  en¬ 
thalpy,  1,  directly  by  starting  from  the  outer  edge  and  proceeding  to 
the  wall. 

As  in  the  explicit  method  the  continuity  equation  is  now  used  to 
determine  the  velocity,  v,  or  V,  across  the  boundary  layer.  Two  methods 
are  used  to  write  the  continuity  equation  as  a  difference  equation.  The 
first  one  is  the  same  as  that  used  in  the  explicit  method  and  is  called 
Method  A.  The  second  one  replaces  the  derivatives  with  difference  quo¬ 
tients  such  that  the  truncation  error  is  of  higher  order  and  is  called 
Method  B.  By  repeated  use  of  the  procedures  outlined  above  and  presented 
in  this  chapter  the  boundary-layer  flow  at  succeeding  steps  downstream 
can  be  determined. 
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In  the  methods  of  solution  of  the  houndaxy- layer  equations  dis¬ 
cussed  above  it  has  been  assumed  that  the  boundary  conditions  are  well 
defined.  Of  course,  at  the  wedl  they  are  well  defined  if  we  assume  no 
slip  as  the  velocity  u  will  be  zero.  The  velocity  v  or  V  is 
usually  taken  as  zero,  except  when  there  is  fluid  injection  and  in  either 
case  no  difficulties  are  introduced.  Also,  either  the  temperature  or 
heat  transfer  distribution  at  the  wall  can  be  easily  handled.  However, 
at  the  outer  edge  the  velocity  and  enthalpy  are  known  from  the  inviscid 
flow  pressure  distribution,  but  there  is  no  definite  location  of  the 
outer  edge.  For  the  boundary- layer  equations  in  the  physical  plane  or 
after  the  Howarth-Dorodnltsyn  transformation  the  velocities  euid  enthalpy 
approach  the  inviscid  flow  values  asymptotically.  There  are  two  methods 
that  can  be  used  in  applying  the  outer  edge  boundary  conditions.  One 
way  is  to  pick  a  line  which  is  parallel  to  the  wall  and  is  sufficiently 
far  away  from  the  wall  that  the  outer  edge  boundary  conditions  are 
valid  there.  Since  the  boimdary  layer  is  normally  thickest  at '.the  down¬ 
stream  limit  of  the  computation,  an  estimate  must  be  made  of  this  thick¬ 
ness.  This  method  is  Inefficient  as  too  many  steps  are  taken  across  the 
boundary  layer  where  the  layer  is  thinner  than  the  maximum  thickness. 

The  second  method,  which  is  used  in  this  paper,  determines  the  edge  of 
the  boundary  layer  by  finding  out  when  the  quantities  determined  across 
the  layer  become  nearly  a  constant.  By  testing  successive  values  of 
the  quantities  calculated  across  the  layer,  it  is  assumed  that  the  edge 
is'  reached  when  the  two  quemtities  are  the  same  to  a  certain  niamber  of 
decimal  places.  The  number  of  decimal  places  required  is  easily  ob¬ 
tained  after  some  experience,  and  this  question  is  illustrated  in 
Chapter  4  with  numerical  solutions.  With  this  method  the  number  of 
steps  across  the  boundary  layer  will  vary  as  the  thickness  varies. 

It  has  been  assumed  that  initial  profiles  of  velocity  u  and  v 
or  Y  and  enthalpy  i  are  known.  Obtaining  these  profiles  is  a  ser¬ 
ious  problem  connected  with  all  the  numerical  schemes  for  solving  the 
boundary- layer  equations.  Since  only  a  local  solution  is  required, 
the  problem  is  considerably  simpler  than  solving  the  boundary- layer 
equations  for  an  arbitrary  body.  Some  of  the  profiles  available  and 
the  onps  used  in  this  paper  are  discussed  later  in  this  chapter  (Section  D). 
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Also  (discussed  in  this  chapter  is  the  stability  of  the  difference 
equations  with  respect  to  round-off  error  and  the  convergence  of  the 
numerical  solution  toward  the  exact  solution.  The  numerical  formxjlas 
used  to  determine  the  boundary- layer  characteristics  of  shearing  stress, 
heat  transfer  and  displacement  thiclmess  sure  also  presented.  Finally, 
the  chapter  is  concluded  with  the  computer  program  used  to  solve  the 
boundary- layer  problem. 

A.  Difference  Quotients 

In  constructing  the  difference  quotients,  the  sketch  below  is  use¬ 
ful  for  reference. 


It  is  assumed  that  the  functions  H(x,y)  and  T(>C;y)  are  known  at  a,  b, 
c,  but  unknown  at  d,  e,  f.  In  the  following  sections  the  difference 
quotients  that  replace  the  partial  derivatives  are  given  as  are  the 
higher-order  terms.  The  higher  order  terms  sure  obtained  using  a  Taylor 
series  expansion  and  are  shown  to  indicate  the  truncation  error  intro¬ 
duced  when  the  first  terms  are  used  to  replace  the  derivatives. 

1.  Explicit.  In  the  explicit  method  the  partial  derivatives 
are  evalmted  at  the  grid  point  "b"  and  the  following  difference  quo¬ 
tients  eo’e  used: 
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All  the  quantities  are  evaluated  at  point  "b"  unless  indicated  other- 
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wise,  and  the  truncation  error  is  of  order  Ax  or  /^y  . 

2.  Implicit . 

a.  Method  I.  This  method  evaluates  the  partial  derivatives 
at  the  grid  point  "e"  since  an  implicit  scheme  is  desired.  The  partial 
dei'ivatives  are  replaced  by  the  difference  quotients  as  indicated  and 
all  quantities  are  evaluated  at  point  "e"  unless  denoted  otherwise. 
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In  the  following,  difference  expressions  for  products  of  derivatives 
are  given.  They  are  chosen  such  that  in  the  products  the  unknown  quan¬ 
tities  appear  linearly  and  their  use  csui  lead  to  linear  difference 
equations . 
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The  truncation  error  in  the  above  difference  quotients  is  of  the  order 
2 

of  Ax  or  Ay  and  is  the  same  as  the  explicit  method. 

b.  Method  II.  This  method  evaluates  the  partial  derivatives 
at  the  grid  point  "k"  since  the  Crank-Nicolson  type  of  difference  equa¬ 
tions  are  desired.  The  partial  derivatives  are  replaced  by  the  follow¬ 
ing  difference  quotients : 
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The  above  difference  quotients  are  chosen  because  the  unknown  quan¬ 
tities  appear  linearly  and  can  result  in  linear  difference  equations. 

All  the  quantities  are  evaluated  at  point  "k"  unless  indicated  otherwise. 

2  2 

The  truncation  error  in  this  case  is  of  order  Ax  and  Ay  and  is  of 
higher  order  in  Ax  than  in  the  previous  cases. 

B.  Difference  Equations 

The  following  sketch  shows  the  designation  of  mesh  points^  in  the 
physical  plane. 


If  the  Howarth-Dorodnitsyn  transformation  is  used,  one  has  a  similar 
grid  in  the  transformed  plane;  but  x  and  y  are  replaced  by  |  and 
■q,  respectively.  It  is  assumed  that  the  velocities,  u  and  v  or  V, 
and  enthalpy,  i,  eire  known  at  the  grid  points  in  the  m^^  column  and  un¬ 
known  at  the  grid  points  in  the  (m+l)^^^  column. 


^This  notation  is  different  from  that  which  hsis  been  used  in  pre¬ 
vious  reports  (Ref.  1  and  11)  on  finite-differqnce  solutions  of  the 
bovindary- layer  equations.  This  change  in  designation  of  mesh  points 
was  instituted  in  order  to  be  consistent  with  the  subscript  notation 
used  in  the  Burroughs  220  Computer  program. 


1.  Explicit.  The  transformed  houndary- layer  equations  are  used 
for  the  explicit  method  since  the  stability  requirements  are  less  severe 
than  when  the  physical- plane  equations  are  employed.  The  boundary- 
layer  equations  (2.29  -  2.31)  become  the  following  when  equations  3>1 
are  used; 
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A-T) 

2S5 


(u. 


m+l,n 


-  u  + 
m.n 


u 


m+l,n-l 


(3.5) 

(3.6) 


The  partial  derivatives  at  (m,n)  have  not  been  actually  replaced  in 
order  to  keep  the  notation  as  simple  as  possible.  However,  they  are 
readily  determined  from  equations  3*1  since  they  do  not  involve  any 
unknowns . 

2.  Implicit.  The  two  implicit  Methods  I  and  II  axe  quite  similar 
and  will  be  discussed  together.  As  indicated  previously,  the  derivatives 
in  the  boundaxy-layer  equations  are  evaluated  at  the  joints  "e"  emd  "k" 
in  the  two  methods  and  correspondingly  the  partial  differential  equations 
must  be  evaluated  at  the  same  place.  Notice  in  equations  2.l8  and  2*30 
that  there  axe  terms  of  the  following  form: 


T 


i 


where  "i"  is  "e"  or  "k",  depending  on  the  method.  If  the 
is  replaced  by  the  difference  quotients  as  given  by  (3* 2a) 
the  above  expression  becomes 


derivative 
or  (3*3a), 
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— f 


Expressions  of  this  type  will  be 
the  difference  equations  will  be 
pansionj  we  can  write  as 


nonlinear  in  the  unknowns  and,  of  course, 
nonlinear  also.  Using  a  Taylor's  ex- 


Ti-i^  +  Ia*  (g)  +  ...  (i-o) 

b 


Therefore,  in  order  to  have  linear  difference  equations,  the  original 
expression  is  written  as 


T 


i 


where  terms  of  order  Ax  have  been  neglected  as  indicated  in  the  above 
expansions . 

Using  the  above  lineatrization  technique  and  the  difference  quo¬ 
tients  (3 -2)  and  (3*3)^  we  can  write  the  momentum  and  energy  equations 
(2.18  -  2.19)  or  (2.30  -  2.31)  of  the  boiindary  layer  as  the  following 
difference  equations: 

(2  <  n  <  N-1  and  m  >  l) 


A-  u,  1+B,  u,  +C,  u,  i+B,  i,  , 
T.  m+l,n-l  1  m+l,n  1  m+l,n+l  1  u+l.n-l 
m,n  '  m,n  '  ’ 


m,n 


m.n 


+  E,  i  ,  F,  1  ,  .  =  G, 

1  m+l,n  1  m+l.n+1  1 

m,n  '  m,n  ’  m,n 


(3.7) 


^2  'm+l,n-l  ^  ®2  ^m+l,n  ^  *^2  '^m+l,n+l  ^  ^2  ^m+l,n-l 

^  win  '  Tnn  '  mn  ' 


m.n 


m.n 


m.n 


®2  ^m+l,n-i  ^2  ^m+l,n+l  °2 
m,n  '  m,n  ’  n 


(3.8) 


where  the  coefficients  vary,  depending  onitiether  Method  I  or  II  is 
being  used.  The  coefficients  for  Method  I,  which  result  from  equations 
(2.18  and  2.19)  with  equations  3.2,  sure  the  following: 
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Method  I  (Physical  Plane) 


m^n 

=  L.  (-pv  +  f  i  - 

(3.9a) 

m,n 

4fL^ 

=  ^  — 

(3.9b) 

m,n 

=  (pv  -  f  i  -  1^) 

1  '  y  ^ 

(3.9c) 

m,n 

=  f  L,  u 

1  y 

(3.9<i) 

m,n 

=  0 

(3.9e) 

m,  n 

=  -  f  Li  u 

^  y 

(3.9?) 

in,n 

2 

=  pu  -  Zac  (p'  +  f  u  i  ) 
e  y  y' 

(3.9g) 

Ag 

m,n 

=  2f  Uy 

(3.10a) 

^2 

m,n 

=  0 

(3.10b) 

^2 

m^n 

=  -  2f  L,  u 

^  y 

(3.10c) 

°2 

m,n 

Pf**  O-P 

Li  (  pv  +  pj,  iy  p^^) 

(3.10d) 

®2 

m^n 

4fL^ 

=  PrAy 

(3.10e) 

F 

2 

m^n 

T  /  ,  2f'  2f  . 

"  1  ^  Pr  y  PrAy 

(3.10f) 

^2 

m,n 

=  pui  +  Ax  (u  p;  -  f  Uy  -  i^) 

(3.10g) 
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where 


Lj^  =  Ax/(2Ay) 

(3.11a) 

i  =  (1  ,  -  i  , )/ (2Ay) 

y  '  m,n+l  m,n-l'''^'' 

(3.11b) 

u  =  (u  1  -  u  . )/ (2Ay) 

y  '  m,n+l  m,n-l'''^^' 

(3.11c) 

In  the  above  relations  all  the  quantities  are  evalmted  at  the  grid  point 
(m,n)  due  to  the  linearization.  Without  the  linearization  eLLl  the  un- 
subscripted  quantities  would  be  evaluated  at  the  grid  point  (m+l^n). 

The  coefficients  for  Method  II,  which  results  from  equations 
(2.30  and  2.31)  with  equations  3*3  are  the  following: 


Method  II  (Transformed  Plane) 


and 


m,n 


m,n 


m.n 


m.n 


m,n 


m.n 


m.n 


^2 


m.n 


m,n 


(3.12a) 

4  FLp 

-in 

(3.12b) 

^2  -  i  -  V 

(3.12c) 

F’  Lo  u,^ 

2  T) 

(3.12d) 

0 

(3.12e) 

'  ^2  N 

(3.12f) 

(3.12g) 

2L2  F  u^ 

(3.13a) 

0 

(3.13b) 
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(3.13c) 


=  -  2  1-2  F  u 
m,n 

4fl 

Eo  =  2u  + 

2m.  n 


^  T  /  ,T  2F  2F'  .  . 

■  2  ^  ■  Prai  -  Pr  1  ) 


=  2u  i  +  Af  (2up'/p  “  V  i  +  77-  i  „) 

m,  n  ^  .  in  Pr 


(3.13<i) 

(3.13e) 

(3.13f) 

(3.13g) 


where 


^2 

=  As/(2  ATI) 

(3.14a) 

^^m,n+l 

i 

m,n- 

i)/(2Atj) 

(3.l4b) 

=  (u  - 

'  m,n+l 

u 

m^n- 

l)/(2An) 

(3.14c) 

^^m,n+l 

2i 

m,n 

+  i  1 

m,n-J 

.)/An^ 

(3.l4d) 

"nr, 

-  - 

2u 

m^n 

+  u  , 

m,n-J 

)lt^ 

(3.l4e) 

In  the  above  relations  all  the  quantities  are  evaluated  at  the  grid  point 
(m,n)  due  to  the  linearization.  Without  the  linearization  all  the  unsub- 
scripted  quantities  would  be  evaluated  at  the  grid  point  (m+l/2,n). 

Two  methods  are  used  to  replace  the  continuity  equation  and  will 
be  called  Methods  A  and  B.  The  following  sketch  of  the  grid  will  be 
useful  in  the  description  of  these  methods  of  writing  the  difference 
equations . 
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Method  A  evaluates  the  derivatives  in  the  continuity  equations  at  point 
"g"  and  uses  the  backward  difference  scheme  for  x-derivatives .  The  dif¬ 
ference  quotients  are 


(pV)y  = 


(pv) 


m+l,n 


(pv) 


m+l,n 


Ay 


,n-l  ^1.2,  . 

^ ^  Ay  (pv) 


yyy 


(pu)^  -  (P^Vn  + 

•  5  (3-15) 

The  continuity  equation  (2.17)  becomes  the  following  when  the  above 
difference  quotients  are  used. 

(Physical  Plane  -  Method  A) 

(3.16) 
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For  Method  B  the  derivatives  are  evaluated  at  point  "h"  as  shown  in  the 
preceding  sketch.  The  difference  quotients  used  in  this  method  are 


k  ("“’xxx  +  S  («)xyjr  + 


(3-17) 


-  (pv) 

m,n  '  'm,n- 


+  -H  Ax^  (pv)  +  Ay^  (pv)  +  .. . 

It  should  be  noticed  that  the  truncation  error  in  these  difference 
quotients  is  of  higher  order  in  the  x-direction  than  in  Method  A.  Us¬ 
ing  the  above  difference  quotients  gives  for  the  continuity  equation 
(2.17)  the  following; 

(Physical  Plane  -  Method  B) 

(3.18) 

^  ^^'^^m+l,n-l  " 

The  continuity  equation  (2.29)  in  transformed  plane  becomes  the 
following  for  the  two  methods  discussed  above. 

(Transformed  Plane  -  Method  A) 


=  V 

m+l,n  m+l^n-1 


An 

2A| 


(^1 


m+l,n  ^m,n  ^  "‘m+l^n-l 


+  u 


u  ,  ) 
m^  n-1' 


(3.19) 
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(Transformed  Plane  -  Method  B) 


V  =  V  -  V  +  V 

m+l,n  m+l,n-l  m,n  m,n-l 


ATI 


(3.20) 


-u  +u,  ,-u  ,) 

A|  m+l,n  m,n  m+l,n-l  m,n-l' 


C.  Method  of  Solution 

In  order  to  solve  the  difference  equations  either  in  the  explicit 
or  implicit  form,  it  is  necessary  to  have  the  houndeiry  conditions  in  a 
form  suitable  for  numerical  computation.  The  velocity  and  enthalpy  at 
the  outer  edge  of  the  boimdary  layer  are  known  in  the  case  where  there 
is  no  Interaction  and  from  (2.23)  become 


u  „  =  u  (x) 
m, N  e '  ' 

(3.21a) 

i  =  i  (x) 
m,N  e'  ' 

(3.21b) 

where  N  is  the  nximber  of  grid  points  across  the  boundary  layer.  At 
the  wall  the  velocity  is  zero  and  the  other  boundary  condition  depends 
upon  whether  the  temperature  of  the  wall  or  the  heat  transfer  is  speci¬ 
fied.  These  boundary  conditions  are  written  as 


or 


V  T  =  0 
m,  1 

(3.22a) 

u  ,  =  0 
m,  1 

(3.22b) 

i  T  =  i  (x) 
m,l  w'  ' 

(3.22c) 

Rs  — —  (i  -  1  ) 

^  Ay  Pr  '  m,2  m,l' 

(  3.22(1) 

The  approximation  for  the  above  derivative  in  (3.22d)  is  rather  crude 
and  this  point  will  be  discussed  later  when  an  insulated  wall  is  con¬ 
sidered.  The  other  condition  that  must  be  given  is  the  velocity  and 


enthalpy  profiles  across  the  boundary  layer  at  the  start.  The  determin¬ 
ation  of  initial  profiles  is  considered  in  Section  D. 

1.  Explicit  Difference  Equations.  The  method  of  solving  the 

explicit  difference  equations  (3.U  -  3-6)  is  apparent.  Since  all  the 
quantities  on  the  right  side  of  the  equations  are  known  from  the  ini¬ 
tial  profiles,  the  velocities  and  enthalpy  are  calculated  directly  at 
each  grid  point  across  the  boundary  layer.  After  starting  at  the  wall 
the  calculations  proceed  until  the  velocity  ^  and  enthalpy 

are  sufficiently  close  to  the  velocity  u^  and  enthalpy  i^  determined 
from  the  external  flow.  Then  the  profiles  across  the  boundary  layer 
are  calculated  at  succeeding  steps  downstream  until  the  desired  distance 
is  obtained. 

2.  Implicit  Difference  Equations.  The  method  of  solving  the  im¬ 
plicit  difference  equations  (3-7)  (3-8)  is  the  same  for  both  impli¬ 

cit  difference  Methods  I  and  II.  Of  course,  the  coefficients  in  the 
difference  equations  will  depend  upon  the  method  being  used.  Obviously, 
the  difference  equations  cannot  be  solved  directly  as  in  the  explicit 
case,  but  one  has  the  problem  of  solving  a  large  number  of  linear  alge¬ 
braic  equations.  The  system  of  algebraic  equations  is  of  a  special 
type  since  a  large  niomber  of  the  coefficients  in  the  complete  system 
are  zero.  Because  of  the  special  form  of  the  equations,  the  following 
relations  exist  (see  Appendix  B,  especially  equations  l4b): 


u  ,  ,  =  yS^}  ,  +  1  u  -  +  ,  1  ,  (3.23a) 

m+l,n-l  m+l,n-l  m+l,n-l  m+l,n  m+l,n-l  m+l,n  '  ' 

i  1  1  1  i-i  n  +  1  i  1  (3.23t>) 

m+l,n-l  m+l,n-l  m+l,n-l  m+l,n  m+l,n-l  m+l,n  ' 

When  the  above  are  substituted  into  the  difference  equations  (3*7)  and 
(3.8),  the  following  is  obtained: 


®1  'V+l,n  ^1  ^m+l,n  ^1  ^1  '^m+l,n+l  ^1  ^m+l,n+l 

m,n  ’  m,n  ’  m,n  m,n  '  m,n 


C3.24a) 
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\+l,n  ^m+l,n  "  ^2 

m,n  ’  '  ’ 


m^n 


m,n  m,n 


^2  %+l,n+l  ■  ^2  ^m+l,n+l 

TYl  n  '  Tn  n  ' 


m,n 


where 


B*  =  +  Aj_ 

m.n  m.n  m.n 


(3-24h) 


*  ”l  4"l,n-l  (3'25>) 

■’  m.n  ’ 


B*  =  B^  +  K 


(2) 


(2) 


^2  -  X.2  -  -2  -mU,n-l  °2  ^m+l,n-l  (3-25t) 

"  “  “  ’  m.n  ’ 


m.n  m.n  m.n 


E*  =  +  A,  -,  +  Di  L  ,1  T 

111  m+l,n-l  1  m+l,n-l 
-  "  --  m  r,  >  m.n  ' 


(3)  (3.25c) 


m.n  m.n  m.n 


^  “2  '■m+ln-1  <3.25d) 

m^n  m,n  m,n  ’  m,n  ’ 


G*  =  G,  -  A,  ,  -  D-  ,  (3.25e) 

111  m+n,n-l  1  m+l,n-l  '  ' 

m,  n  m,  n  ’  m,  n  ’ 


G*  =  G2  -  Ag  K 
m,n  m,n 


(1) 


(1) 


2  “  "^2  "2  *m+l,n-l  '  ^2  ^m+l,n-l  (3*25f) 

“  “  ’  m.n  ' 


Solving  equation  (3.24a)  and  (3.24b)  for  and  ^m+l  n 


u  =  + 

m+1, n  m+1, n 


u  + 

m+1, n  m+1, n+1 


K 


(3)  i 

■m+l,n  m+l,n+l 


(3.26a) 


1  ,  .l(\>  * 

m+l,n  m+l^n 


u  ,  + 

m+l,n  m+l,n+l 


.(3)  i 
"'m+l,n  m+l,n+l 


(3.26b) 


where 


4ii,n  =  -  “I 


(3.27a) 


m+l,n 


=  A  (E*  C2 


C,  E*) 


1  2^m,n 
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(3.27b) 


-(3) 


K:i,n  =  ^  ^2  - 


(3.27c) 


(1) 


^  (^1  ^Pra.n 


(3.27d) 


(2) 


r,  =  ^  (Cl  B§  -  B*  C„)^  ^ 

nrt-l^n  1  2  1  2'in,  n 


(3.27e) 


(3) 


^  (^1  ^2)m,n 


(3.27f) 


where 


A  =  1/(BI  E.  -  EJ 


After  we  apply  the  boundary  condition  at  the  wall  (3.22b),  equa¬ 
tions  (15B)  in  Appendix  B  require  the  following: 


=  K<^>  .  k(7  = 

m,  1  m,  1  m,  1 


(3.28) 


Similarly,  the  boundary  condition  (3.22c) 


and  equations  (15B)  give 


Specified  wall  temperature  (3.29) 


while  boundary  condition  (3.22d)  gives 


m,  1 
m,  1 
m,  1 


Ay 

Pr  Sf 


0 

1 


Specified  heat  transfer  (3.30) 


.Now  the  necessary  relations  exist  to  solve  the  implicit  finite-difference 


equations,  but  the  method  needs  to  be  clarified. 

.(1)  j^(2)  ^(3)  ^(1)  ^(2) 


The  quantities  K'  , 
m,n’ 


,  L'“', 

m^n'  m,n'  m,n'  ^,n' 


and  L 


(3) 

m,n 


can  be  de¬ 


termined  across  the  boundary  layer  by  using  the  following  procedure: 


(a)  Perform  the  following  steps  at  the  first  point  away  from 
the  wall: 

(1)  Calculate  C^,  E^,  and  from 

equations  (3.9a-  -3.9g)  or  (3.12a  -  3.12g). 

(2)  Calciilate  C^,  ^2 

equations  (3.10a  -  3.10g)  or  (3.13a  -  3.13g). 

(3)  Using  the  results  from  the  previous  steps  and  the 

boundary  conditions  (3.28)  and  (3.29)  or  (3.30), 
calculate  B^*,  E^*,  Gj_*,  and  G2*  as 

given  by  equations  (3.25a  -  3-25^). 

(k)  Using  previous  results,  calculate 

and  with  equations 

(3.27a  -  3.27f). 

m 

(b)  Now  the  procedure  outlined  in  step  (a)  can  be  repeated  at  the 

second  point  away  from  the  wall  using  results  obtained  at  the 

/ 

first  point.  The  above  procedxire  is  repeated  until  the  boun¬ 
dary  layer  is  traversed  and  all  values  of 
and  are  determined. 


Now  the  velocity,  u,  and  enthalpy,  i,  are  determined  across  the 
boundary  layer,  starting  at  the  outer  edge.  Knowing  the  K' s  and  L's 
across  the  boundary  layer,  we  can  use  equations  (3.26a)  and  (3.26b)  to 
solve  for  the  velocity  and  enthalpy  utilizing  the  boundary  conditions 
(3.21a)  and  (3.21b)  at  the  outer  edge. 

As  yet  the  continuity  equation  has  not  been  required  in  order  to 
calculate  the  velocity,  v  or  V,  due  to  the  linearization  of  the  dif¬ 
ference  equations.  Before  the  computations  can  proceed  downstream, 
the  velocity,  v  or  V,  must  be  determined  across  the  boundary  layer. 
Starting  at  the  wall  calculate  the  velocities  v  or  V  from  equation 
(3.16)  or  (3.18)  for  the  physical  plane  and  equation  (3.1^)  ot  (3.20) 
for  the  transformed  plane,  respectively.  Now  the  velocity  and  enthalpy 
profiles  can  be  determined  at  succeeding  steps  downstream. 
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D.  Initial  Profiles 


In  most  of  the  problems  solved  using  numerical  schemes^  initial  pro¬ 
files  are  obtained  from  similarity  solutions  of  the  boundary- layer  equa¬ 
tions.  However,  Wu  (Ref.  42)  has  proposed  the  following  type  of  initial 
profiles  for  leading  edges  or  stagnation  flow: 

(1)  The  velocity  u  and  enthalpy  i  are  freestream  values  at  all 
the  grid  points  across  the  layer  except  at  the  wall. 

(2)  At  the  wall  the  velocity  u  is  zero  and  the  enthalpy  corres¬ 
ponds  to  the  wall  temperature. 

(3)  The  velocity  V  is  assxmied  zero  at  all  the  grid  points  across 
the  boundary  layer. 

Since  this  method  would  be  very  advantageous  for  starting  numericeO. 
computations,  it  has  been  investigated  further.^ 

To  study  the  possibility  of  using  the  Wu  type  initial  profiles,  the 
boundary- layer  flow  from  the  leading  edge  of  a  flat  plate  at  Mach  number 
9.6  was  solved  using  the  explicit  method.  A  linear  viscosity  law  was 
used  so  that  the  results  could  be  compared  to  similarity  solutions  which 
were  obtained  from  Low  (Ref.  26).  The  displacement  thickness  of  the 
boundary  layer  has  been  used  to  show  the  influence  of  using  Wu  type  ini¬ 
tial  profiles  on  the  niimerical  solutions.  The  displacement  thickness 
along  a  flat  plate  for  two  sizes  of  mesh  normal  to  the  wall,  AH,  and 
two  sizes  of  mesh  parallel  to  the  wall,  A§,  is  shown  in  Fig.  la.  It 
seems  that  the  numerical  solutions  can  be  close  to  the  exact  solution 
at  a  sufficient  distance  downstream  if  the  proper  step-sizes  are  chosen. 
However,  simply  reducing  the  step-sizes  does  not  always  improve  the  re¬ 
sults  . 

In  the  Wu  type  initial  profiles  it  has  been  assumed  that  V  is 
zero  at  the  leading  edge,  but  the  similar  solution  gives  V  equal  to 
infinity.  From  physical  consideration  the  assumption  V  =  0  is  reason¬ 
able;  but  from  the  mathematical  point  of  view,  the  similar  solution 

^It  should  be  mentioned  that  such  profiles  are  not  correct  for  stag¬ 
nation  flow,  e.g.,  as  the  mesh  size  approaches  zero.  For  incompressible 
flow  at  a  stagnation  point  the  displacement  thickness  is  finite,  while 
the  Wu  profiles  wo\iLd  give  zero  for  this  quantity. 
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Figure  1,  Plat  Plate  Plow  with  Wu  Type  Initial  Profiles  at  the 

Leading  Edge 
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is  correct.  In  other  words,  the  horn dary- layer  equations  are  not  valid 
near  the  leading  edge,  which  is  well  known.  Therefore,  in  the  next 
test  the  velocity,  V,  for  the  initial  profile  was  assumed  at  all  grid 
points  across  the  boundary  layer  eqml  to  the  value  given  by  the  simi¬ 
larity  solutions  at  the  outer  edge  and  at  the  first  step  from  the  lead¬ 
ing  edge.  At  the  wall  the  velocity  V  was  assumed  zero  as  \iBual.  The 
influence  on  the  displacement  thickness  by  assuming  V  finite  at  the 
leading  edge  is  shown  in  Fig.  lb.  In  this  case  the  nimierlcal  solution 
agrees  closely  with  the  similar  solution. 

From  the  above  discussion  it  is  concluded  that  Wu  type  initial 
profiles  can  give  reasonable  results  if  the  proper  mesh  sizes  are  choqen. 
Even  better  results  are  obtained  when  the  velocity  V  is  given  the 
appropriate  value  rather  than  zero.  Since  it  is  difficult  to  ascertain 
the  required  mesh  sizes  and  velocity  V,  the  Wu  type  initial  profiles 
must  be  used  with  extreme  care. 

To  obtain  initisLl  profiles  it  seems  better  to  use  similai'  solutions 
rather  than  the  Wu  type  of  initial  profiles.  There  are  three  types  of 
initial  bovindary- layer  profiles  which  should  be  sufficient  for  most 
practical  problems.  The  three  types  are  flat  plate,  wedge,  eind  stagna¬ 
tion  flow.  For  a  complete  discussion  of  similar  solutions  of  the  boun¬ 
dary-layer  equations,  see  Chapter  8  of  Hayes  and  Probstein  (Ref.  l4). 

For  supersonic  speeds  with  an  attached  shock  wave,  the  flat  plate  and 
wedge  profiles  are  the  same  (i.e.,  constant-pressure  profiles).  For 
the  case  of  Prandtl  niimber  0.72  and  the  linear  viscosity  law,  the  ini¬ 
tial  profiles  are  easily  obtained  from  Low  (Ref.  26) .  When  the  Sutherland 
viscosity  law  is  used,  the  initial  profiles  can  be  obtained  using  Crocco's 
method  of  solution  as  applied  by  Van  Driest  (Ref.  4o). 

For  two-dimensional  stagnation- point  flow  the  profiles  that  are 
available  are  based  on  the  assumption  that  Prandtl  n\jmber  is  one  or  a 
linear  viscosity  law  is  used.  Since  the  boundary- layer  eqviations  for 
the  stagnation  point  reduce  to  two  simultaneous  nonlinear  ordinary  dif¬ 
ferential  equations,  it  is  possible  to  solve  these  equations  numerically. 
Because  these  equations  are  of  a  "two-point  boundary- value  problem" 
type,  they  are  tedious  to  solve,  but  the  Sutherland  viscosity  law  and 
Pr  =  0.7  can  be  \ised. 
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In  the  verification  of  the  implicit  finite-difference  scheme, 
initial  profiles  for  the  flow  along  a  flat  plate  are  used.  A  linear 
viscosity  law  is  assumed  since  the  profiles  are  easily  obtained  from 
Low  for  this  case.  In  order  to  obtain  the  profiles  in  the  notation  of 
this  paper,  certain  relations  must  be  established  with  Low's  results. 
Where  the  derivation  of  these  relations  become  involved,  the  procedure 
is  presented  in  Appendix  A.  The  necessary  relations  to  obtain  the 
initial  profiles  in  the  physical  plane  are  the  following: 


where 


y  =  Q  g(\o^) 

^  =  i  (\ow) 

^  =  ^e  e’  (\ow) 
u  Q 


®'(Vov)  '  1  =(\ov) 


K  = 


2.0748 


w 


-  1)  -  ^  (0.8477) 


Q  =  2 


Pe  % 


(3.31a) 

(3.31t)) 

(3.31c) 

(3.31<i) 

(3.32a) 

(3.32b) 

(3.32c) 

(3.32d) 


The  relations  used  to  obtain  the  initial  profiles  in  the  transformed 
plane  are  the  following: 


n  -  Pg  Q  ^2iow 

(3.33a) 

‘■'tW 

(3.33t) 

i  =  le  S'  ( 

(3.33c) 
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(3-33d) 


V  >  tv,  f  (-.!„)  -  f(v„)l 

Some  of  the  initial  profiles  are  given  in  Chapter  4  along  vith  the 
numerical  problems  solved. 

E.  Stability  and  Convergence  of  the  Difference  Schemes 

All  of  the  difference  schemes  presented  in  this  paper  are  consistent 
or  are  a  formal  approximation  to  the  partial  differential  equations .  A 
scheme  is  consistent  if  the  difference  between  the  partial  differential 
equation  and  difference  equation  goes  to  zero  as  the  mesh  size  approaches 
zero.  In  other  words,  the  truncation  error  goes  to  zero  as  the  step- 
sizes  go  to  zero,  which  is  easily  seen  from  the  difference  quotients 
used  in  the  various  schemes .  The  faster  the  truncation  error  goes  to 
zero,  the  more  accurate  is  the  difference  approximation  and  more  likely 
is  the  numerical  solution  to  be  a  good  approximation  to  the  exact  solu¬ 
tion.  However,  consistency  does  not  imply  that  the  solution  of  the  dif¬ 
ference  equations  tends  to  the  solution  of  the  partial  differential  equa¬ 
tions  as  the  step-sizes  go  to  zero.  For  this  to  be  true  the  difference 
equations  must  be  convergent.  For  equations  of  the  boundary-layer  type, 
there  is  no  completely  satisfactory  mathematical  analysis  to  show  the 
convergences  of  the  difference  schemes.  Fliigge-Lotz  and  Yu  (Ref.  ll) 
have  made  the  most  significant  investigation  of  the  convergence  of  the 
explicit  difference  scheme.  The  convergence  of  the  difference  equa¬ 
tions  can  be  studied  numerically  by  varying  the  mesh  sizes.  This  type 
of  verification  of  the  convergence  of  the  difference  schemes  is  pre¬ 
sented  in  Chapter  along  with  the  numerical  examples  solved. 

Since  there  are  round-off  errors  in  the  computations,  the  numeri¬ 
cal  error  between  the  exact  and  numericeuL  solutions  of  the  difference 
equations  must  be  investigated.  This  problem  is  referred  to  as  the 
stability  of  the  numerical  scheme.  By  making  several  assumptions  and 
using  the  von  Neumann  method  of  stability  analysis,  we  can  obtain  an 
approximate  estimation  of  step^size  for  stability.  Previous  experience 
(Ref.  1)  indicates  that  the  momentum  equation  dominates  in  determining 
the  stability  of  the  boundary- layer  equations.  Therefore,  only  the 
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momentum  equation  will  be  considered  for  the  stability  of  the  boundary- 
layer  equations.  If  the  pressure  gradient  term  is  assumed  zero  and  a 
linear  viscosity  law  is  used,  the  momentum  equation  (2.30)  in  the  trans 
formed  plane  becomes 


uu^  +  V  u^ 


F  u, 


(3*31+) 


Next  it  is  assumed  that  the  mesh  is  sufficiently  small  so  that  the  co¬ 
efficients  in  the  above  equation  vary  slightly  and  are  considered  con¬ 
stants.  The  above  equation  then  can  be  written  as 


uu^  +  V  u^ 


(3.35) 


where  the  bars  indicate  the  quantities  are  to  be  considered  constants. 
If  the  derivatives  are  replaced  with  the  explicit  difference  quotients 
(3-1)  and  the  von  Neumann  analysis  is  applied,  the  usual  stability  re¬ 
quirement  is  obtained 


Stability  Parameter  =  -  —  <  1  (3*36) 

u  An 

The  bars  have  been  omitted  because  the  local  values  of  the  quantities 
will  be  used.  Since  the  velocity  is  smallest  near  the  wall,  the  first 
grid  point  away  from  the  wall  will  give  the  largest  value  for  the  left- 
hand  side  of  the  inequality  (3*36) • 

The  flow  along  a  flat  p].ate  at  Mach  number  9*6  has  been  solved 
using  the  explicit  method  with  a  grid  size  large  enough  to  cause  in¬ 
stability.  The  calculated  displacement  thickness  is  given  in  Fig.  2 
where  the  instability  occurs  at  approximately  |  =  O.05.  Also  shown  in 
this  figure  is  the  stability  parameter  (3*36)  at  the  first  grid  point 
away  from  the  wall.  The  value  of  the  stability  parameter  is  approxi¬ 
mately  two  when  the  Instability  occurs,  rather  than  one  as  relation 
(3.36)  requires.  Besides  the  assumptions  mentioned  previously,  some 
of  this  difference  can  be  attributed  to  the  fact  that  the  known  boimdary 
conditions  at  the  wall  improve  stability  (see  Ref.  I5). 
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Figure  2.  Displacement  Thickness  for  Unstable  E:5)licit  Difference 
Scheme  Solution  vith  Variation  of  Stability  Parameter 
Shovn 
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Since  the  simple  analysis  above  gives  a  reasonable  estimate  of 
stability  requirements,  the  same  type  of  analysis  is  now  applied  to  the 
second  implicit  difference  scheme.  Substituting  the  difference  quo¬ 
tients  (3-3)  into  the  momentum  equation  (3*35)  gives  the  following  dif¬ 
ference  equation: 


a,  u 


where 


+  b,  u 


and 


Ct  U  T  T  +  d-l 

1  m+l,n+l  1 

'^m,n-l  '''  ®1  ^m,n^^l  '^n,n+l  ^ 

(3.36) 

a^  =  dj^  =  -r  -s 

(3.37a) 

b^  =  1  +  2s 

(3.37b) 

<=1  ■  A  ■ 

(3.37c) 

e^  =  -1  +  2s 

(3.37(i) 

r 

4u  An 

(3.38a) 

2u  An 

(3.38b) 

Applying  the  von  Neumann  method  of  stability  analysis,  we  substitute 

ox  iPy  OinAx  ip(n-l)Z^  .  .  .  .  /_  ^  ... 

u  =  e  e  ''=e  e  '  into  equation  (3»3d),  and  letting 

I  =  e  ,  gives 


I  = 


+  (d^+f^)  cos  PAy  -  i(d^-f^)  sin  PAy 
bi  +  (a^^+c^)  cos  PAy  -  i(a^-c^)  sin  PAy 


(3-39) 


Using  relations  (3-37).>  we  can  write  the  above  as 


The  above  can  be  written  in  the  following  form: 

(1  -  0)(1  +  0)  -  -  129 

where 

0  =  2s(l  -  cos  PAy) 

9  =  2r  sin  PAy 

The  condition  for  stability  is 


111  <  1 

l|i"<l 


From  equation  (3-^l)  the  following  is  obtained; 

I  |2  (l.0f  (l+0f  +  26^  (1+/)  + 

(1+0)^  (1+0)^  +  29^  (l+0f  + 


If  we  consider  the  case  when  u  >  0,  then  0  >  0,  If  0  =  0,  then 
(3.44)  becomes 

for  all  e.  If  0  >  0,  then  the  following  relations  exist; 

(l+0f  >  (l-0f 

(1+0)^  >  1+/ 

Using  these  relations  with  (3*44),  it  is  easy  to  see  that 


(3.41) 

(3.42a) 

(3.42b) 

(3.43a) 

(3.43t) 

(3.44) 

equation 

(3.45a) 

(3.45t) 


• — ! 


for  all  9  when  0  >  0.  Therefore,  the  implicit  Method  II  with  u  >  0  is 
^stable  regardless  of  the  mesh  size  as  shown  in  the  above  analysis.  The 
numerical  examples  investigated  in  Chapter  ]y  corroborate  this  result. 

F.  Formulas  for  Shear  Stress,  Heat  Transfer,  and  Displacement  Thickness 

The  formula  for  the  shear  stress  in  the  physical  plane  at  the  wall 
has  been  presented  in  equation  (2.45),  'but  now  it  is  necessary  to  ex¬ 
press  the  derivative  in  this  expression  numerically.  The  velocity  u 
near  the  wall  may  be  expressed  with  sufficient  accuracy  as 

u  =  7i  y  +  7^  y^  +  7^  (3-46) 

Taicing  the  derivative  of  the  above  with  respect  to  y  and  setting 
y  =  0,  we  can  write  the  fonnula  for  shear  stress  at  the  wall  (2.45) 

\  ‘  ’■l  (3.1.7) 

The  value  of  7^  can  be  determined  by  evaluating  equation  (3.46)  at  the 
three  grid  points  nearest  the  wall  and  solving  the  three  resulting  equa¬ 
tions  for  7^.  When  the  result  for  7^  is  substituted  into  equation  (3.^7)^ 
the  following  is  obtained  for  the  shear  stress  at  the  wall  in  the  physi¬ 
cal  plane  s 

^  =  ^  ^^2  -  9  ^^3  +  2  u^)  (3.48a) 

An  expression  for  the  shear  stress  at  the  wall  in  the  transformed  plane 
is  obtained  in  a  similar  manner  from  equation  (2.46)  and  is 

^  ^2  -  9  ^3  ^  ^4)  ^3. 48b) 

The  formula  for  the  heat  transfer  at  the  wall  has  been  presented  in 
equation  (2.49)  an<i  again  it  is  necessary  to  replace  the  derivative. 

If  the  same  procedure  is  followed  as  above,  except  the  velocity  is 
replaced  by  the  enthalpy  in  all  the  expressions,  the  following  formula 
is  obtained  for  evaluation  of  the  heat  transfer  at  the  wall  in  the 
physical  plane; 
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(3.49a) 


%  "  6Ay  Pr  i-2  ^  ^3  "  ^  h'^ 

The  corresponding  expression  for  the  transformed  plane  is 

'5v  -  1^  -  18  +  9  1,  -  2  I4)  (3.1.gb) 


In  order  to  determine  the  displacement  thickness  as  given  by  equation 


(2„53)  a  trapezoidal  rule  of  integration  is  used.  Since  u  is  zero  at 
the  wall,  the  dimensionless  displacement  thickness  in  the  physiceil  plane 
can  be  written  as 


6*  =  Ay 


(3.50a) 


Using  the  trapezoidal  rule  again,  we  write  the  displacement  thickness 
in  the  transformed  plane  as  given  by  equation  (2.54)  as 


6*  =  An 


e 


(3.50b) 


To  determine  the  displacement  thickness,  a  better  integration  method  is 
Simpson's  rule  which  gives  the  following  relation  in  the  transformed 
plane : 


G.  Computer  Program 

Due  to  the  lengthy  and  repetitious'  nature  of  the  implicit  finite- 
difference  schemes  for  solving  the  boundary- layer  equations,  the  problem 
was  programmed  for  the  IBM  65O  and  Burroughs  220  digital  computers.  The 
flow  diagram  for  the  basic  computer  program  for  the  explicit  difference 
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scheme  without  interaction  is  given  in  Fig.  3a.  For  the  implicit  finite 

difference  schemes  the  explicit  subroutine  for  u  ^  and  i  ,  ,  as 

m+l,n  m+l,n^ 

indicated  in  Figure  3a,  is  replaced  by  the  subroutine  for  computing 
%+l  n  ^m+1  n  "the  implicit  scheme  (see  Figure  3b).  The  flow 

diagram  can  be  divided  into  many  subroutines,  as  shown  in  the  diagram; 
however,  some  operations  are  not  as  clearly  separated  as  indicated. 

Given  below  is  the  computation  procedure  which  will  relate  previous 
formulas  with  the  flow  diagram. 

1.  Program  Constants.  This  consists  of  computing  various  parameters 
which  are  used  frequently  but  do  not  change  =  Also  computed  are 
u^,Po^,  i^,  and  T^  from  equations  (2.36  -  2.39)* 

2.  Compute  p  and  p'  .  The  pressure  and  pressure  gradient  are 

_ ^m+1 _ ^m+1 

computed  at  the  points  m+1.  The  formulas  or  values  for  these  quan¬ 
tities  are  given  \inder  the  specific  problem  being  solved. 

3-  Compute  Exterior  Flow  Quantities .  The  exterior  flow  quantities 
Ug,  ig,  pg,  Tg,  and  are  computed  using  formulas  (2.4o  -  2.44). 

4.  T  <0.  This  is  a  check  to  be  sure  that  there  is  no  laminar 
w 

separation. 

5*  Compute  u  ,  and  i  ^  .  The  velocity  and  enthalpy  across  the 

_ m+l,n _ m+l,n 

boundary  layer  at  mesh  points  m+1  are  calculated.  For  the  method 
of  solution  and  equations  used,  see  Chapter  III,  Section, Cf2. 

6 .  Test  for  Edge  of  Boundary  Layer  While  Computing  u  and  i .  The 

following  test  that  is  described  is  applicable  only  when  the  veloc¬ 
ity  and  enthalpy  of  the  exterior  flow  at  any  x  are  independent  of  y 
In  the  implicit  calculation  proced\ire  it  is  necessary  to  know  how 
far  to  calculate  the  K^^^'s  and  L^^^'s  (i  =  1,2,3)  across  the  boun¬ 
dary  layer.  The  typical  variation  of  and  across  the 

boundary  layer  is  shown  in  Fig.  4.  At  the  outer  edge  of  the  boun¬ 
dary  layer  these  quantities  become  independent  of  the  distance 
away  from  the  wall.  Therefore,  to  stop  the  calculation  of  the 
K^^^'s  and  L^^^'s  (i  =  1,2,3)  after  the  boundary  layer  has  been 
transversed,  two  consecutive  values  of  are  compared  to  see  if 

the  difference  between  them  is  less  than  some  small  quantity  e. 

The  value  of  e  is  determined  by  the  desired  accuracy  of  the 
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Figure  3a.  Flow  Diagram  for  Explicit  Finite-Difference 
Schemes  Without  Interaction 
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Figure  k.  Typical  Variation  of  and  Across  the 

Boundary  Layer 


-  5*t  - 


calculations  and  is  discussed  further  in  Chapter  IVwith  the  numerical 
examples.  If  the  difference  between  the  K^^^'s  is  less  than  e,  the 
difference  between  two  consecutive  L^^^'s  is  compared  with  e.  When 
this  difference  is  also  less  than  e,  the  computations  proceed  to 
the  next  step.  This  procedure  applies  for  computations  in  both  the 
physical  and  transformed  planes. 

Add  u^  and  1^  at  Edge  of  Boundary  Layer.  This  portion  of  the  pro¬ 
gram  adds  u^  and  i^  at  the  mesh  points  beyond  the  last  points  com¬ 
puted.  This  step  is  required  so  that  the  number  of  mesh  points 
across  the  boundary  layer  will  not  decrease  by  one  at  each  step 
downstream. 

Compute  or  The  velocity  v  or  V  is  computed 

across  the  boundary  layer  using  equations  (3.I6  or  3.I8)  and  equa¬ 
tions  (3*19  or  3*20),  respectively. 

Test  for  Edge  of  Boundary  Layer  While  Computing  v  or  V.  This 
time  the  test  for  the  edge  of  the  boundary  layer  is  simply  to  have 
the  same  number  of  mesh  points  across  the  boundary  layer  as  for  the 
u  and  i  computations. 

Compute  6*  and  Wall  Quantities.  The  displacement  thickness,  6*, 
is  computed  from  the  appropriate  equation  of  (3.50a)  to  (3.50c). 

The  wall  quantities  and  q^  are  computed  from  equations  (3.48) 
and  (3*49),  respectively.  The  shear  stress  parameter  C^y Re  ' 
and  heat  transfer  parameter  StyRe^  are  computed  using  equations 

(2.47)  and  (2.50). 


CHAPTER  IV 


ILLUSTRATIONS  OF  THE  IMPLICIT  DIFFERENCE 
SCHEMES  WITHOUT  INTERACTION 


The  results  for  a  variety  of  boundary-layer  problems,  using  the 
implicit  difference  schemes  developed  in  the  previous  chapter,  will  now 
be  presented.  Although  Implicit  Method  II  is  superior,  it  has  the  dis¬ 
advantage  of  using  the  transformed  plane.  Since  Implicit  Method  I 
solves  the  boundary-layer  equation  in  the  physical  plane,  it  is  inter¬ 
esting  to  consider  examples  solved  by  this  method  and  understand  its 
features.  Four  problems  are  investigated  using  Method  I,  and  two 
problems  with  Method  II.  In  addition,  the  flow  near  the  leading  edge 
of  a  flat  plate  is  solved  by  both  methods. 

A.  Implicit  Method  I 

Five  problems  were  solved  with  an  IBM  65O  computer  using  Implicit 
Method  I,  which  solves  the  boundary-layer  equations  in  the  physical 
plane.  The  first  two  problems  involve  the  flow  along  a  flat  plate  where 
the  initial  profiles  are  at  one  reference  length  from  the  leading  edge 
(x^  =  1.00).  In  example  one  the  wall  temperature  is  assumed  constant, 
but  with  a  value  such  that  the  wall  is  heated.  The  second  flat  plate 
problem  assumes  the  wall  is  insulated  (zero  heat  transfer).  In  both 
of  these  examples  a  linear  viscosity  law  is  used  so  that  a  direct  compari¬ 
son  can  be  made  with  similar  solutions  of  the  boundary-layer  equations. 
Then  in  example  three  the  flow  along  a  wall  with  a  ramp  pressure  gradient 
is  investigated  and  is  compared  with  numerical  results  of  Baxter.  Since 
Baxter  uses  Sutherland's  viscosity  law  in  this  example,  the  present 
example  uses  the  same  viscosity  law.  To  study  the  influence  of  the 
viscosity  law  on  the  boundary-layer  solution,  the  fourth  example  con¬ 
siders  the  flow  along  a  flat  plate  with  the  same  conditions  as  example 
one  except  Sutherland's  viscosity  law  is  used.  Finally,  the  flow  near 
the  leading  edge  of  a  flat  plate  (x  -•  O.Ol)  is  Investigated,  since  this 
case  will  be  considered  with  displacement  thickness  interaction  in 
Chapter  V. 
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1.  Flat  Plate  FIov  With  Constant  Wall  Teniperature .  The  boundary- 
layer  flow  along  a  flat  plate  whose  surface  is  heated  and  temperature 
is  constant  has  been  investigated  by  Flugge-Lotz  and  Yu  (Ref.  ll).  An 
explicit  finite-difference  scheme  was  used  by  Flugge-Lotz  and  Yu  for 
solving  the  boundary-layer  equations  in  the  physical  plane.  For  this 
example  the  numerical  solution  became  unstable  after  a  few  steps  from 
the  initial  profiles.  Also,  the  solution  to  this  problem  can  be 
obtained  from  similarity  results  and  is  readily  obtained  from  Low 
(Ref.  26).  Since  the  similar  solution  is  based  upon  a  linear  viscosity 
law  (Eq.  2.21),  the  same  assumption  is  used  in  the  numerical  solution. 

The  following  fluid  properties,  exterior  flow,  and  wall  quantity, 
are  used  for  this  example: 


M  =  3  (U.la) 

00  « 

T*  =  389.99°^  (^-It) 

CD 

7  =  l.k  (^.Ic) 

Pr  =  0.72  (^.Id) 

S*  =  216°  R  (^.le) 


The  pressure  distribution  has  not  been  discussed  previously.  For  a 
flat  plate  without  Interaction  the  following  relations  are  used  : 


Pe  = 


K 


=  0 


(4.2a) 

(4.2b) 


As  a  consequence  of  the  above  relations,  the  exterior  flow  quantities 
at  the  edge  of  the  boundary  layer  are  the  same  as  the  corresponding 
quantities  at  the  free  stream  conditions.  The  initial  profiles  are 
obtained  from  equations  (3.31®'’  -  3«31'i)  and  are  shown  in  Figure  5* 

This  problem  has  been  solved  with  Implicit  Method  I  and  no  stability 
difficulties  are  encountered.  In  these  calculations  step  -sizes  up  to 
one  hundred  times  greater  than  the  size  employed  by  Yu  have  been  used. 
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The  velocity  profiles  u  and  v  and  the  enthalpy  profile  i  obtained 
from  the  numerical  solution,  after  ten  steps  downstream  of  the  initial 
profiles,  are  shown  in  Figure  5*  The  similar  solution  to  the  boundary- 
layer  equations  at  the  same  distance  downstream  is  also  given. 

In  the  numerical  computations  the  parameter  e  and  the  step-size 
must  be  chosen.  In  Chapter  III,  Section  0.6  the  computer  program  test 
for  the  edge  of  the  boundary  layer  was  discussed  and  certain  quantities 
had  to  be  less  than  e.  The  value  of  e  was  determined  by  decreasing 
its  magnitude  until  there  was  no  influence  on  the  displacement  thick¬ 
ness.  For  this  problem  €  =  0.00001.  Several  step-sizes  have  been 
used  to  solve  this  problem.  The  effect  of  changing  Ay  and  Ax  on 
the  boundary-layer  characteristics  is  shown  in  Figures  6  and  7  respec¬ 
tively.  These  characteristics  were  obtained  from  equations  (3»^a), 
(3.49a),  and  (3.50a).  The  boundary- layer  characteristics,  as  given  by 
similar  solution,  are  also  shown  in  these  figures.  For  the  range  of 
step-sizes  investigated  there  is  very  little  influence  on  the  conver¬ 
gence  of  the  numerical  solution  to  the  exact  solution.  It  should  be 
noticed  that  the  scales  in  these  figures  are  greatly  expanded  and  the 
differences  between  the  numerical  and  similar  solutions  are  actually 
very  small.  Since  the  initial  profiles  are  only  accurate  to  four 
decimal  places,  there  are  initially  errors  in  the  boundary-layer  cal¬ 
culations  and  characteristics.  These  errors  appear  to  decrease  as  the 
computations  proceed, 

2.  Flat  Plate  Flow  With  Zero  Heat  Transfer.  At  the  wall  either 
the  wall  temperature  or  the  heat  transfer  can  be  specified.  Since  the 
previous  example  considered  the  specified  wall  temperature,  this  example 
takes  the  case  when  the  heat  transfer  is  zero  (insulated  wall).  This 
example  is  the  same  as  the  problem  in  the  previous  section  except  for 
the  wall  boundary  condition.  The  fluid  properties  and  exterior  flow 
are  given  by  expressions  (4.1a  -  4.1e),  and  the  relations  (4.2a  -  4.2b) 
are  also  applicable  for  the  pressure  distribution.  The  initial  profiles 
are  shown  in  Figure  9  and  were  obtained  from  equations  (3.31)* 

The  wall  boundary  condition  for  the  case  of  a  specified  heat 
transfer  is  given  in  equation  (3.22d).  The  values  for  the  coefficients 
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due  to  this  houndary  condition. 

ni,l  111,1'  ni,l 

(3*30).  For  the  case  of  zero  heat  transfer,  this 
gives  the  following  results: 


are  given  in  equation 
boundary  condition 


m,  1  m, 2 

(4.3a) 

L<\)  =  0 
m,  1 

(4.3’o) 

l'")  =  0 

m,  1 

(4.3c) 

=  1 

m,  1 

(4.3<i) 

As  previously  mentioned,  the  derivative  in  the  boundary  condition  (3.22d) 
has  been  replaced  by  a  rather  crude  approximation.  By  using  two  mesh 
points  beyond  the  wall  in  writing  the  difference  quotient,  a  better 
approximation  for  the  derivative  is  obtained.  For  an  insulated  wall 
the  following  relation  results: 


i  =  —  i  -  —  i 

^m,l  3  ni,2  3  m,3 


(^.4) 


If  n  is  set  equal  to  two  and  u  ^  is  set  equal  to  zero  (boundary 

condition  3* 22b)  in  equations  3*7  and  3*8^  then  u  _  can  be  eliminated 

m,  3 

from  the  two  resulting  equations.  As  a  result  of  these  operations,  the 
following  equation  is  obtained  : 


(B1C2  -  %,2  ^  ^°l‘^2  ■  ®2^1^m,2  ^m,l  ^®1^2  "  ^2^1  ^ 2  ^m,2 


+  ^m,3  ^^1^2  ■  ^2^1  ^m, 2 


(4.5) 


Substituting  equation  (4,4)  into  the  above  equation  (4,5)  and 
rearranging  gives : 


[3(FiC2  -  F^C^)  -  (D^C^  - 

+  (B^C2  -  %2  f^V2  • 


■  ^1^2 ^m, 2 

“  ^2^1 ^^m, 2  ^m,2 


(4.6) 
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By  comparing  coefficients  in  the  above  equation  with  equation  (3.26b), 
the  following  is  obtained: 

^  <®1=2  - 

^  [(^2  -  ^=1)  *  ‘‘('■A 

where 

l/A,  .  [3(F^C2  -  F^C,)  -  (D,02  -  D20^)1„^2 

The  flow  along  an  insulated  flat  plate  has  been  solved  using  the 
different  boundary  conditions  (4.3)  and  (4.7).  In  Figure  8  the  wall 
enthalpy  is  shown  when  the  two  different  boundary  conditions  are  used, 
and  the  similarity  result  is  also  shown.  The  profiles  across  the 
boundary  layer  at  x  =  1.10  are  shown  in  Figure  9*  The  numerical  solu¬ 
tion  with  boundary  conditions  (4.7)  are  compared  with  the  similar  solu¬ 
tion  in  this  figure.  The  boundary-layer  characteristics  for  the 
insulated  flat  plate  are  presented  in  Figure  10.  The  numerical  results 
with  the  two  different  boundary  conditions  are  compared  with  the  similar 
solutions  in  this  figure.  From  these  results  it  is  apparent  that 
boundary  conditions  4.3  are  indeed  crude.  Therefore,  to  approximate 
the  enthalpy  gradient  at  the  wall  it  is  necessary  to  use  at  least  three 
grid  points  as  was  done  for  the  boundary  conditions  (4.7). 

3.  Ramp  Pressure  Gradient.  In  order  to  investigate  the  influence 
of  the  boundary- layer  pressure  gradient  term  on  the  numerical  solution, 
the  flow  along  a  wall  with  a  constant  wall  temperature  and  a  ramp 
pressure  gradient  is  considered.  This  problem  has  been  solved  by 
Baxter  and  Fliigge-Lotz  (Ref.  l),  using  an  explicit  finite-difference 
scheme  with  the  Crocco  transformed  boundary-layer  equations.  The  fluid 
properties,  exterior  flow,  and  wall  quantities  are  the  same  as  the  first 
example  and  are  given  in  relations  (4.1).  Since  Baxter  uses  Sutherland's 
viscosity  law  (Eq.  2.22),  it  is  also  used  in  this  example.  The  pressure 
and  pressure  gradient  are  the  following 
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(4.7a) 

(4.7b) 

■  ^2®1^L,2  (4.7c) 


V, 


Boundary  Conditions  4.7) 


Sinin^  Solution 


Figure  10.  Insulated  Flat  Plate  Boundary-Layer  Characteristics 


"  Poo  +  2  >  1) 

=  5  (x  -  1)  p^u^  (x  >  1)  (4.8b) 


for  the  ramp  pressure  gradient  case.  Flow  along  a  flat  plate  is  assumed 
upstream  of  x  =  1,  and  the  initial  profiles  given  in  Figure  5  are  used 
to  start  the  calculations.  Rather  than  use  equation  (2.4o),  the  small- 
disturbance  form  of  this  equation  was  used.  The  simpler  equation  is; 


00 


(4.9) 


In  the  finite -difference  solution  Ay  =  1.0  and  Ax  =  0.001  or 
0.0004.  The  latter  value  of  the  step-size  was  used  by  Baxter  in  order 
to  have  a  stable  numerical  solution.  The  parameter  e  used  in  testing 

-4 

for  the  edge  of  the  boundary  wa^'  increased  to  10  for  this  example. 

The  boundary-layer  characteristics  for  a  ramp  pressure  gradient 
are  given  in  Figure  11.  The  shearing  stress  parameter  was  obtained  from 
equations  (2.47)  and  (3.48a),  while  the  heat  transfer  parameter  was 
obtained  from  equatiohs  (2.50)  and  (3*49a).  The  displacement  thickness 
parameter  is  the  ratio  of  the  dimensionless  displacement  thickness  at 
any  x  to  the  value  at  x  =  1.0.  The  displacement  thickness  is  obtained 
from  formula  (3.50a).  The  step-size  must  be  the  same  size  (Ax  =  0.0004) 
as  used  by  Baxter  in  order  to  have  reasonable  agreement  with  his  results. 
Since  the  boundary-layer  profiles  are  changing  very  rapidly,  the  trun¬ 
cation  error  is  more  important  than  stability  considerations  in  deter¬ 
mining  the  step-size.  As  the  truncation  error  in  Baxter's  and  the 
present  method  are  of  the  same  order,  it  is  not  surprising  that  the 
required  step-sizes  are  approximately  the  same  for  this  example. 

It  should  be  remembered  that  the  Initial  profiles  are  based  upon  a 
linear  viscosity  law,  while  the  numerical  results  use  Sutherland's  vis¬ 
cosity  law.  Near  the  start  of  the  computation  of  the  boundary-layer 
characteristics,  the  influence  of  changing  the  viscosity  law  can  be 
seen  in  Figure  11. 
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4,  Flat  Plate  Flow  With  Sutherland's  Viscosity  Lav.  In  example  3 
the  initial  profiles  were  based  upon  a  linear  viscosity  law,  while  the 
numerical  computation  used  Sutherland's  viscosity  law.  In  order  to 
investigate  the  influence  of  such  a  procedure  on  the  boundary-layer 
characteristics,  example  1  is  solved  again  using  Sutherland's  viscosity 
law. 

In  Figure  12  the  variation  of  viscosity  with  enthalpy  for  the  two 
viscosity  laws  is  illustrated.  For  a  heated  wall  and  large  Mach  number 
the  linear  viscosity  law  becomes  a  poor  approximation  for  the  more  exact 
Sutherland's  law. 

The  boundary-layer  characteristics  for  this  example  are  given  in 
Figure  13.  The  results  for  the  similar  solution  with  a  linear  viscosity 
law  are  from  Low  (Ref.  26),  while  those  with  Sutherland's  viscosity  law 
are  from  Van  Driest  (Ref.  40).  Two  step-sizes  were  used  in  the  numerical 
computations  and  these  results  are  presented  in  this  figure.  The 
boundary-layer  characteristics  which  are  obtained  from  the  numerical 
solution  with  Sutherland's  viscosity  law  tend  to  approach  the  similar 
solution  from  Van  Driest.  The  computations  with  the  small  step- size 
seem  to  be  approaching  the  results  of  Van  Driest  faster  than  the  other 
numerical  solution.  As  this  type  of  problem  is  considered  again  in 
Section  B  example  1  of  this  chapter  (See  Figure  20 ),  the  computations 
have  not  been  extended  further  downstream. 

5.  Flow  Near  the  Leading  Edge  of  a  Flat  Plate.  This  example 
considers  the  flow  near  the  leading  edge  of  a  flat  plate  when  the  free- 
stream  Mach  number  is  9 •6.  The  linear  viscosity  law  is  used  so  that 
the  results  can  be  compared  to  the  similar  solution  easily  obtained  from 
Low  (Ref.  26).  The  following  fluid  properties,  exterior  flow,  and  wall 
quantity,  are  assumed  : 


M  =  9.6 

00 

(4.10a) 

T*  =  82.34°R 

00 

(4.10b) 

7  =  1.4 

(4.10c) 

Pr  =  0.72 

(4.10d) 
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Figure  12.  Vlecoslty  Lavs 


S*=  198.^  R 


(i+.lOe) 

(i+.lOf) 


T  /T  ,=  1 
w'  ad 

The  wall  temperature  is  specified  by  equation  (4.10f)  such  that  there 

is  no  heat  transfer  in  this  example.  The  pressure  distribution  is  still 

given  by  equations  (4.2)  since  this  is  flow  along  a  flat  plate.  The 

initial  profiles  are  obtained  from  equations  (3.31a  -  3-31<i)  and  are 

shown  in  Figure  l4.  In  this  problem  the  test  for  the  edge  of  the 

-I4. 

boundary  layer  uses  e  =  10  . 

In  this  problem  the  initial  profiles  are  taken  at  x^  =  0.01  and 
the  calculations  proceed  to  x  =  0.10.  This  is  equivalent  to  starting 
at  x^  =  1.00  as  done  in  example  1  (Section  A.l)  but  continuing  the 
computations  until  x  =  10.00.  Therefore  a  relatively  large  step-size 
is  required  if  the  finite -difference  scheme  is  going  to  be  of  any 
practical  value  in  this  example. 

Two  step-sizes  have  been  used  to  solve  this  problem  and  the 
boundary- layer  displacement  thickness  is  used  to  illustrate  the  resulting 
solutions  in  Figure  15 .  The  similar  solution  is  shown  in  this  figure 
also.  This  figure  shows  that  the  grid  size  is  too  large  and  hence  the 
truncation  error  has  significant  influence  on  the  results.  As  the  grid 
size  Ax  is  reduced,  the  solution  appears  to  be  converging  to  the 
similar  solution.  When  the  numerical  solution  is  to  be  sufficiently 
accurate,  the  required  grid  size  would  become  prohibitively  small. 

Another  factor  has  been  Introduced  which  requires  the  grid  size  to 
be  small  in  this  problem.  In  order  to  approximate  the  profiles  as  shown 
in  Figure  l4  with  sufficient  accuracy,  it  is  necessary  to  have  a  small 
grid  size  Ay  near  the  outer  edge  of  the  boundary  layer.  Shown  in 
this  figure  is  the  numerical  and  similar  solution  at  one  step  downstream 
from  the  initial  profiles.  The  numerical  solution  is  greatly  different 
from  the  similar  solution  near  the  outer  edge  of  the  boundary  layer. 

The  problem  of  reducing  the  truncation  error  and  of  increasing  the 
step-size  near  the  outer  edge  of  the  boundary  layer  requires  going  to 
Implicit  Method  II. 
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B.  Implicit  Method  II 


Three  examples  were  solved  with  a  Burroughs  220  computer  using 
Implicit  Method  II,  which  solves  the  boundary-layer  equations  in  the 
transformed  plane.  The  first  example  is  the  same  as  the  last  case 
considered  in  the  previous  section.  This  example  illustrates  that 
Method  II  has  a  much  smaller  truncation  error  than  Method  I.  Other 
effects  studied  in  this  example  are:  the  two  methods  of  replacing  the 
continuity  equation,  integration  formula  for  the  displacement  thickness, 
value  of  €  for  test  at  edge  of  boundary  layer,  iteration  of  the  solu¬ 
tion  at  each  step  downstream,  and  type  of  viscosity  law.  In  the  second 
example  the  flow  downstream  of  a  transpiration  cooled  region  is  investi¬ 
gated.  The  wall  in  the  region  of  interest  has  “  0.5.  This 

problem  and  the  insulated  wall  case  have  been  solved  by  Howe  (Ref.  17) 
using  the  numerical  method  of  Baxter  and  Flugge-Lotz  (Ref.  l).  The 
insulated  wall  case  has  been  solved  by  Pallone  (Ref,  3^)  with  another 
numerical  method.  This  example  allows  an  indirect  comparison  between 
these  numerical  methods  for  solving  the  boundary-layer  equations.  The 
third  example  Illustrates  that  the  implicit  method  can  be  used  to  solve 
the  boundary-layer  flow  along  a  wall  with  a  variable-  temperature.  The 
specific  case  considered  is  a  -wall  with  a  "hot  spot"  as  this  problem 
has  been  solved  by  Baxter  and  Flugge-Lotz  (Ref,  l). 

1.  Flow  Near  the  Leading  Edge  of  a  Flat  Plate,  The  same  problem 
as  given  in  Section  A. 5  of  this  chapter  is  now  investigated.  The  dis¬ 
placement  thickness  of  the  boundary  layer  is  used  to  Illustrate  the 
influences  of  several  parameters  in  Figure  l6.  In  all  parts  of  the 
figure  the  similar  solution  of  the  boundary-layer  equations  is  given  to 
indicate  the  desired  result.  In  part  (a)  the  result  using  Implicit 
Method  I  for  this  problem  as  given  in  Section  A. 5  is  presented  for 
convenience.  Also  in  this  part  of  the  figure  the  displacement  thickness 
is  shown  for  the  case  when  the  problem  is  solved  in  the  physical  plane 
but  with  the  difference  quotients  used  in  Implicit  Method  II.  By  using 
the  new  difference  quotients,  the  truncation  error  has  been  greatly 
reduced.  However,  when  the  profiles  are  examined  for  this  solution, 
the  numerical  result  still  has  some  error  near  the  outer  edge. 
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8lBil«r  flolutloa  (law) 


In  Figure  1613  the  influence  of  transforming  the  boundary- layer 
equations  is  illustrated.  The  displacement  thickness^  as  a  result  of 
solving  these  transformed  equations,  does  not  seem  to  he  Improved. 

This  numerical  solution  does  approximate  the  similar  solution  profiles 
closer  than  the  physical  plane  solution,  especially  near  the  outer  edge. 
It  should  he  remembered  that  all  solutions  have  used  Method  A  for  the 
continuity  equation.  In  Figure  l6c  the  Influence  of  using  the  two 
methods  of  replacing  the  continuity  equation  is  shown.  When  Method  B 
is  used  with  Implicit  Method  II  the  best  numerical  solution  of  the  dis¬ 
placement  thickness  is  obtained.  For  this  case  the  nximerlcal  profiles 
are  a  good  approximation  for  the  similar  solution  results. 

Several  other  items  that  have  a  smaller  Influence  on  the  solution 
of  this  problem,  have  been  investigated.  In  order  to  illustrate  the 
Influence  of  these  items,  the  displacement  thickness  is  divided  by  the 
similar  solution  displacement  thickness  at  the  corresponding  value  of 
Xo  This  displacement  thickness  ratio  is  given  in  Figure  17  for  the 
various  cases  studied.  The  numerical  values  have  been  connected  by 
straight  lines  to  clarify  the  figure.  The  first  item  investigated  was 
the  influence  of  the  integration  formula  used  in  evaluating  the  dis¬ 
placement  thickness.  The  displacement  thicknesses  computed,  using  the 
trapezoidal  rule  (Eq.  and  Simpson's  rule  (Eq.  3-50c)>  give 

essentially  the  same  result  for  this  problem.  Since  Simpson's  rule  is 
a  better  formula  in  general  and  is  only  slightly  more  complicated  than 
the  trapezoidal  rule,  Simpson's  rule  is  used  with  Implicit  Method  II. 

The  influence  of  €  on  the  displacement  thickness  ratio  is  shown  in 
Figure  17.  From  this  figure  it  is  impossible  to  ascertain  the  better 
value  of  €  to  use.  To  study  the  effect  of  e  further,  the  velocity 
profile  near  the  outer  edge  is  given  in  Figure  l8.  This  figure  gives 
the  velocity  profile  at  x  =  0,125  for  the  numerical  solution  with 
two  values  of  €  and  the  similar  solution  is  also  given.  This  figure 
shows  that  the  smaller  value  of  €  results  in  a  velocity  profile  which 
is  a  better  approximation  to  the  similar  solution. 

Returning  to  Figure  17,  another  item  that  has  been  investigated  is 
the  effect  of  iterating  Implicit  Method  II  at  each  downstream  step. 

In  order  to  obtain  linear  difference  equations  certain  quantities  at 
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(m  +  -g-)  were  replaced  by  known  quantities  at  m..  (See  Chapter  III, 
Section  B.2),  After  the  boundary-layer  equations  are  solved  as  usual 
for  the  unknown  quantities  at  (tn  +  l),  the  equations  are  resolved  with 
the  quantities  at  (m  +  replaced  by  the  average  of  the  corresponding 
quantities  at  (m)  and  (m  +  l).  By  iterating  in  such  a  manner,  the 
truncation  error  should  be  reduced.  This  is  the  case  as  illustrated,  in 
Figure  17  where  the  solution  was  iterated  once  at  each  step.  Since  this 
iteration  almost  doubles  the  computation  time,  it  would  be  equivalent 
to  halving  the  step-size  downstream.  The  results  for  the  displacement 
thickness  ratio  for  this  smaller  step-size  are  given  in  Figure  17. 
Although  this  result  is  closer  to  the  similar  solution  than  previous 
cases  considered,  it  is  not  as  close  as  the  iterated  solution.  However, 
the  iterated  method  has  a  small  oscillation  that  does  not  occur  when 
the  usual  method  is  used. 

The  comparison  between  Implicit  Method  II  (Transformed  Plane) 
and  the  similar  solution  is  now  completed  by  considering  profiles  and 
shearing  stress  at  the  wall.  For  this  comparison  the  results  are 
obtained  by  using  continuity  equation  Method  B  with  e  =  0,001, 

Ax  =  0,0025,  and  Ar\  =  0.00023288.^  The  initial  profiles  in  the  trans¬ 
formed  plane  are  given  in  Figure  19  and  are  obtained  from  equations 
(3* 33)0  The  numerical  and  similar  profiles  at  x  =  0.125  are  shown 
in  this  figure.  The  two  velocity  and  enthalpy  profiles  obtained  by 
both  methods  agree,  as  expected  from  the  displacement  thickness  results. 
In  Figure  20  the  shearing  stress  at  the  wall  from  the  numerical  solution 
is  given.  This  was  obtained  from  eqxxation  (3.^b).  The  shearing  stress 
from  the  similar  solution  with  linear  viscosity  (Ref.  26)  is  presented 
in  this  figure  for  comparison  with  the  numerical  solution.  Since  the 
initial  profiles  for  the  numerical  computation  were  determined  with  a 
linear  viscosity  law,  the  same  law  is  normally  used  in  the  numerical 
computations.  Here,  however,  the  numerical  solution  has  been  performed 
with  Sutherland's  viscosity  law  and  the  result  for  the  shearing  stress 
is  shown  in  Figure  20a.  The  skin-friction  parameter  for  this  example 


■^his  number  for  Atj  occurs  as  Atilo.j,  is  taken  as  a  convenient 
value  in  determining  the  initial  profiles  from  equations  (3-33)- 
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Figure  20.  Shearing  Stress  and  Shin  Friction  for  the  Flov  Near  the 
Leading  Edge  of  a  Flat  Plate  (M  =  9.6} 


is  given  in  Figure  20b.  The  numerical  solutions  obtained  with  the 

linear  and  Sutherland's  viscosity  law  are  shown.  In  addition  the  skin- 

friction  parameter,  as  obtained  from  similar  solutions  for  the  two 

viscosity  laws,  is  given.  The  similar  solution  result  with  a  linear 

viscosity  ■^law  is  obtained  from  Low  (Ref.  26),  while  the  result  with 

Sutherland's  law  is  obtained  from  Yoimg  and  Janssen  (Ref.  43).  The 

latter  result  is  not  obtained  using  the  same  conditions  as  the  numerical 

solution.  The  skin-friction  parameter  obtained  from  Young  and  Janssen 

was  determined  with  the  following  different  conditions:  T*  =  100°R, 

M  =  9«57»  and  variable  Prandtl  number.  (See  equations  4.10  for  condi- 
00 

tlons  in  present  problem. )  For  the  case  with  a  linear  viscosity  law, 
the  larger  error  near  the  start  can  be  attributed  to  the  relatively 
large  grid  size  in  this  region.  When  Sutherland's  viscosity  law  is 
employed,  the  skin-friction  parameter  near  the  start  is  over  corrected 
as  a  result  of  the  error  introduced  by  the  initial  profiles.  The  skin- 
friction  parameter  then  becomes  nearly  constant  but  at  a  slightly  higher 
value  than  the  similar  solution.  Some  of  this  difference  can  be 
accounted  for  by  the  different  conditions  employed  in  the  two  problems. 
The  heat  transfer  for  this  case  is  not  presented  since  it  is  approxi¬ 
mately  zero  as  the  wall  temperature  was  specified  to  correspond  to  an 
Insulated  wall. 

2.  Flow  Downstream  of  a  Transpiration  Cooled  Region.  This 
example  is  concerned  with  the  flow  downstream  of  a  transpiration  cooled 
region.  The  flow  conditions,  fluid  properties,  and  wall  condition 
are  the  following: 


M  =  3 

00 

(4.11a) 

T*  =  389^99°  R 

00 

(4.11b) 

y  —  1.4 

(4.11c) 

Pr  =  0.72 

(4. lid) 

S*  =  216°  R 

(4. lie) 

'ad  = 

(4. Ilf) 
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The  wall  for  this  case  Is  taken  as  a  flat  plate;  hence,  the  pressure 
gradient  is  zero.  The  injection  region  extends  from  the  leading  edge 
to  X  =  1.0  and  is  such  that  the  boundary- layer  flow  is  described  by 
the  similar  solutions  of  Low  (Ref.  27).  Hence,  the  injection  velocity 
at  the  wall,  v^,  varies  parabolically  from  the  leading  edge.  The 
initial  profiles  are  obtained  from  equations  (3.33).  The  f  and  g 
functions  are  evaluated  with  f(0)  =  -  » —  (p  v  /p  u  )  ]! p  u  x/p  =  -1.0 

from  Low  (Ref.  27).  The  initial  profiles  are  based  upon  a  linear  vis¬ 
cosity  law  and  are  presented  in  Figure  21,  At  the  wall  the  transformed 
normal  velocity  V  is  not  zero  for  the  Initial  profile  as  shown  in 
Figure  21c.  Since  the  wall  is  solid  downstream  of  x  =  1.0,  the  velocity 

V  is  zero  for  x  >  1.0.  In  the  test  for  the  edge  of  the  boundary  layer 

in  the  numerical  computations,  e  =  10  .  The  viscosity  law  used  in 

these  computations  is  indicated  when  the  results  are  presented. 

In  Figure  21  the  boiuidary- layer  profiles  for  a  linear  viscosity 
law  are  presented  at  x  =  2,00  where  the  numerical  computations  were 
terminated.  There  is  very  little  change  in  the  thickness  of  the 
boundary  layer  in  the  transformed  plane.  However,  there  has  been  an 
appreciable  change  in  the  shape  of  the  profiles.  The  quantity  pv  is 
presented  in  Figure  21c  as  it  has  a  more  physical  significance  while 

V  does  not . 

The  skin-friction  ratio  has  been  employed  to  illustrate  the  solu¬ 
tion  of  the  boundary-layer  equation  for  various  conditions  in  this 
problem.  The  skin-friction  ratio  is  obtained  by  dividing  the  local 
skin-friction  coefficient  by  the  following  flat-plate  value  based  upon 
a  linear  viscosity  law: 

(C„)  =  0.6641  l/c/Re  '  (4.12) 

^  f(0)  =0  ^ 

The  skin-friction  ratio  obtained  from  Implicit  Method  II  with  several 
step-sizes  of  Ax  is  presented  in  Figure  22a.  As  the  step-size  Ax 
is  decreased,  the  truncation  error  should  be  reduced  and  the  numerical 
solution  should  be  a  closer  approximation  to  the  exact  solution  (see 
Ref.  32).  When  Ax  =  0.004,  the  results  Indicate  that  the  truncation 
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Figure  21.  Boundary-layer  Profiles  for  Flow  Downstream  of 

Transpiration  Cooled  Region  -  Lineau*  Viscosity  Law 
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Figure  22.  Boundary-Layer  Characteristics  for  Flow 

Downstream  of  Transpiration  Cooled  Region 


error  attributed  to  Ax  is  reasonably  small.  The  step-size  across 
the  boundary  layer  has  been  varied  also.  For  the  two  values  of  A  T) 
there  is  an  appreciable  change  in  the  skin-friction  ratio.  However, 
the  truncation  error  due  to  A  is  small  for  the  smaller  step- size 
across  the  boundary  layer.  It  has  been  estimated  that  the  result  for 
the  case  when  Ax  =  0.01  and  A  t)  =  0.055192  will  change  about  1  per¬ 
cent  near  x  =  2.0  as  the  grid  size  goes  to  zero.  In  Figure  22a  the 
numerical  solution  by  Howe,  which  uses  the  Fliigge-Lotz  and  Baxter  method 
as  given  in  Ref.  10,  is  presented.  This  result  is  based  upon  Suther¬ 
land's  viscosity  law  and  the  same  conditions  (4.11)  are  used  in  this 
example.  These  two  methods  give  results  within  5  percent,  which  is 
very  good  in  comparison  to  other  approximate  methods  for  solving  this 
problem  as  discussed  by  Howe  (Ref.  17)-  However,  one  would  expect  the 
results  to  be  even  closer  as  the  methods  should  be  close  to  the  exact 
solution  if  the  step-size  is  sufficiently  small. 

Since  the  initial  profiles  in  this  problem  are  determined  using  a 
linear  viscosity  law,  this  problem  has  been  solved  with  the  same  law 
to  eliminate  any  initial  error.  The  result  for  the  skin-friction  ratio 
is  given  in  Figure  22b  when  Implicit  Method  II  and  Howe's  method  are 
employed.  The  results  for  two  grid  sizes  with  Howe's  numerical  method 
are  presented  in  this  figure.  As  the  grid  size  is  reduced,  the  skin- 
friction  ratio  is  decreased  -vrtien  Howe's  method  is  used  while  the  skin- 
friction  ratio  is  Increased  when  the  implicit  method  is  used.  Since  the 
results  in  Figure  22b  are  reasonably  close,  the  two  methods  would 
become  very  close  as  the  grid  size  is  decreased.  Since  there  is  a  dis¬ 
continuity  in  certain  quantities  at  the  Initial  profiles,  there  are 
unusually  large  errors  introduced  near  the  start  of  the  computations. 

The  effect  of  these  errors  can,  perhaps,  contribute  to  any  remaining 
difference  between  the  two  methods  with  a  very  small  grid  size. 

When  the  viscosity  law  is  changed  there  is  no  significant  variation 
in  the  skin-friction  ratio  as  obtained  by  Howe's  method.  The  implicit 
method  Indicates  approximately  a  5  percent  decrease  in  the  skin-friction 
ratio  when  Sutherland's  viscosity  law  is  employed.  Because  of  this 
discrepancy,  this  problem  has  been  solved  using  the  explicit  method  as 
described  in  Chapter  III  Section  B.l.  The  results  of  the  explicit  method 
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are  In  agreement  with  the  results  obtained  with  the  implicit  method. 

It  seems  that  the  method  employed  by  Howe  Is  insensitive  to  the  viscosity 
law  employed.  This  peculiarity  of  the  numerical  solution  of  the  Crocco 
form  of  the  boundary- layer  equations  is  discussed  further  in  the  next 
example . 

In  Figure  22c  the  heat  transfer  parameter  for  this  problem  is 
presented.  It  was  obtained  by  dividing  the  nondimens ional  local  heat 
transfer  rate  for  this  problem  by  the  same  quantity  for  a  flat  plate 
with  a  constant  wall  temperature.  The  latter  quantity  called 
(q)f(0)_o  is  the  following  for  this  problem  and  is  based  upon  a  linear 
viscosity  law: 


^ -  (‘*•"3) 

8 . 05^58 

In  Figure  22c  the  heat  transfer  parameter  results  obtained  with  the 
linear  and  Sutherland's  viscosity  law  with  Implicit  Method  II  are 
presented.  The  result  for  this  parameter  as  obtained  by  Howe  is  also 
shown  and  is  computed  with  Sutherland's  viscosity  law.  However,  his 
result  is  in  closer  agreement  with  the  implicit  result  obtained  with 
the  linear  viscosity. 

The  displacement  thickness  is  presented  in  Figure  22d  for  the 
Implicit  results  when  the  two  viscosity  laws  are  employed.  There  is 
first  a  decrease  in  displacement  thickness  after  the  transpiration 
cooled  region.  This  kind  of  behavior  has  been  predicted  by  Pallone 
(Ref.  3^)  for  this  type  of  problem  with  an  insulated  wall. 

As  discussed  previously,  these  examples  provide  an  indirect 
method  of  comparing  several  of  the  numerical  schemes.  However,  this 
comparison  can  be  only  approximate  as  different  digital  computers  were 
used,  computation  time  is  very  sensitive  to  the  accuracy  obtained,  and 
slightly  different  problems  were  solved.  A  problem  similar  to  the 
present  example,  except  the  wall  is  insulated,  has  been  solved  with  a 
linear  viscosity  law  by  Pallone  (Ref.  3^)»  The  same  problem  has  been 
solved  by  Howe  (Ref.  17)  with  Sutherland's  viscosity  law.  The  two 
methods  give  approximately  the  same  result  for  the  skin-friction  ratio. 
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Although  one  would  expect  some  differences  between  the  two  results 
because  of  the  different  viscosity  laws,  the  peculiarity  of  the  Howe 
method  discussed  earlier  can  account  for  the  results  being  nearly  the 
same.  The  computation  of  Pallone  was  performed  with  an  IBM  70^  computer 
in  20  minutes.  Since  the  computations  went  from  x=1.0  to  x=4.0 
and  the  IBM  704  is  approximately  twice  as  fast  as  the  Burroughs  220,  the 
Pallone  method  would  probably  require  about  20  minutes  with  the  Burroughs 
220  for  the  present  problem.  The  computation  of  Howe  was  performed 
with  an  IBM  65O.  computer  a&d  required  7  l/2  hours.  For  the  present 
problem  on  the  Burroughs  220  the  Howe  method  would  require  about  90 
minutes.  When  Implicit  Method  II  is  used,  from  11  to  22  minutes  are 
required  for  the  results  obtained  for  this  problem  with  the  Burroughs 
220  computer.  When  Sutherland's  viscosity  law  is  employed,  approxi¬ 
mately  50  percent  more  computation  time  is  required  than  the  results 
obtained  with  a  linear  law.  The  more  accurate  results  with  the  smaller 
grid  size  take  the  larger  computation  time  indicated.  When  the  same 
grid  size  is  used  as  in  the  implicit  scheme,  the  explicit  method  requires 
about  27  percent  less  computation  time.  However,  the  truncation  error 
appears  to  be  approximately  the  same  when  the  grid  size  Ax  for  the 
explicit  scheme  is  2.5  times  smaller.  Therefore,  the  Implicit  Method  II 
requires  less  computer  time  than  the  eijqplicit  scheme  with  the  same 
accuracy.  From  the  above  considerations,  this  problem  can  be  solved 
with  Implicit  Method  II  as  efficiently  as  any  of  the  other  methods. 

3.  Flat  Plate  Flow  with  "Hot  Spot".  To  illustrate  that  the 
implicit  scheme  can  be  used  to  solve  boundary-layer  flows  with  a 
variable  wall  temperature,  the  problem  solved  by  Baxter  and  Flugge-Lotz 
(Ref,  1)  of  flow  along'a  flat  plate  with  a  "hot  spot"  is  now  considered. 
The  same  conditions  are  used  so  that  a  comparison  of  results  can  be 
made.  These  conditions  are: 


M  =  0.5 

00 

(4.l4a) 

T*  =  389 -99°^ 

00 

(4.l4b) 

7  =  1.4 

(4,l4c) 

Pr  =  0.72 
S*  =  2160r 


(4.l4d) 

(4.l4e) 


The  following  equation  is  used  to  specify  the  wall  enthalpy  for 

1  <  X  <  1.0256: 


i  = 
w 


2  +  0.04 


3/  \3 

X  -  l\  / 1.0256  -  xl 


Vo. 0128/  V  0.0128 


ad 


(4.15a) 


For  X  >  1.0256,  the  wall  enthalpy  is 

7-1 


i  =  2 

w 


1  + 


0.845 


i  =  2i  , 
e  ad 


(4.15b) 


As  this  is  a  flat  plate  problem,  the  pressure  gradient  is  zero.  When 
the  linear  viscosity  law  is  used,  the  constant  C  is  determined  using 
the  wall  temperature  upstream  of  the  "hot  spot". 

The  boundary-layer  characteristics  for  this  problem  are  presented 
in  Figure  23.  The  skin-friction  and  heat-transfer  parameters  have  been 
determined  with  Implicit  Method  II  with  the  two  viscosity  laws.  The 
result  of  Baxter  and  Flugge-Lotz  is  presented  in  this  figure  and  was 
computed  using  Sutherland's  viscosity  law.  In  Figure  23a  the  value  of 
the  skin-friction  parameter  for  a  flat  plate  without  a  "hot  spot"  is 
presented.  This  result  with  a  linear  viscosity  law  was  obtained  from 
Low  (Ref.  26),  while  the  result  with  Sutherland's  viscosity  law  was 
obtained  from  Crocco  (Ref.  9)»  Far  downstream  from  the  "hot  spot"  the 
boundary-layer  characteristics  should  approach  the  flat  plate  similar 
solution  results.  The  pecularity  of  the  Crocco  finite-difference  scheme 
of  numerical  solution  appears  in  this  problem  again.  The  skin-friction 
parameter  obtained  by  Baxter  and  Flugge-Lotz,  which  uses  the  Crocco 
method,  does  not  appear  to  approach  the  similar  solution  of  Crocco. 
However,  the  implicit  solution  with  Sutherland's  viscosity  law  appears 
to  be  approaching  the  Crocco  similar  solution.  When  the  linear  viscosity 
law  is  employed  with  the  implicit  method,  the  variation  of  the  skin- 
friction  parameter  is  similar  to  the  Baxter  and  Flugge-Lotz  result.  The 
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shift  in  the  skin-friction  occurs  because  the  step-size  A  t|  across  the 
layer  is  not  sufficiently  small. 

In  Figure  23b  the  implicit  results  for  the  heat-transfer  parameter 
are  compared  with  the  Baxter  and  Flligge-Lotz  results.  The  viscosity 
law  employed  in  the  implicit  solution  has  only  a  small  effect  on  the 
results.  The  implicit  method  predicts  a  higher  peak  with  a  shift  down¬ 
stream  when  compared  with  the  Baxter  and  Flugge-Lotz  result.  This 
difference  is  insignificant  when  compared  to  other  methods  of  predicting 
the  heat  transfer,  and  the  conclusions  of  Baxter  and  Flligge-Lotz  are 
still  valid. 

In  order  to  investigate  the  closeness  of  the  numerical  and  Crocco 
similar  solutions  further  downstream,  the  flow  along  the  same  flat  plate 
without  the  "hot  spot"  has  been  solved.  The  skin-friction  parameter  is 
presented  in  Figure  23c.  Two  step-sizes  have  been  employed  and  the  two 
viscosity  Laws  have  been  used  with  Implicit  Method  II.  For  the  skin- 
friction  parameter  with  the  linear  viscosity  law,  reasonable  agreement 
is  obtained  with  the  similar  solution  when  the  step-size  is  sufficiently 
small.  When  Sutherland's  viscosity  law  is  employed,  the  implicit  scheme 
indicates  a  slightly  larger  value  of  the  skin-friction  than  the  similar 
solution.  The  present  problem  has  been  solved  with  Sutherland's  vis¬ 
cosity  law  by  Baxter  and  Flugge-Lotz  (Ref.  l),  but  the  computations 
only  extend  a  short  distance  downstream.  This  computation  has  been 
extended  further  downstream,  and  the  result  for  the  skin-friction  is 
presented  in  Figure  23c. 

The  solution  is  coming  closer  to  the  Crocco  similar  solution,  but 
the  process  is  very  slow.  Therefore,  the  Baxter  and  Flligge-Lotz  method 
appears  to  be  insensitive  to  the  viscosity  law  employed  in  the  examples 
Investigated. 

For  the  Implicit  Method  II  solution  in  Figure  23c  with  the  larger 
value  of  Ax,  an  oscillation  of  the  skin-friction  parameter  occurs. 

This  is  attributed  to  an  error  being  introduced  by  the  linearization 
of  the  difference  equations  and  propagated  when  the  continuity  equation 
is  replaced  by  Method  B  (see  equation  3*20).  After  the  two  methods  of 
replacing  the  continuity  equation  have  been  studied,  the  following 
statements  can  be  made:  (l)  For  Method  A  errors  in  V  are  only 
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propagated  indirectly  by  the  influence  on  u  and  i  at  the  next  step. 
(2)  For  Method  B  an  isolated  error  in  V  will  oscillate  in  sign  as 
the  computations  proceed  and  the  indirect  Influence  will  also  occur. 
Therefore,  Method  B  can  give  results  with  a  small  oscillation,  but  with 
less  over-all  error. 
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CHAPTER  V 


NUMERICAL  SOLUTION  OF  THE  BOUNDARY-LAYER  EQUATIONS 
WITH  DISPLACEMENT  THICKNESS  INTERACTION 

In  this  chapter  only  the  boundary-layer  flow  along  a  flat  plate 
with  a  sharp  leading  edge  will  be  Investigated.  However,  the  displace¬ 
ment  thickness  interaction  between  the  boundary-layer  and  external  flow 
will  be  included  in  the  solution.  The  reasons  for  neglecting  other 
parameters  that  contribute  to  the  interaction  between  the  boundary-layer 
and  external  flow  were  discussed  in  Chapter  I. 

In  Hayes  and  Probsteln  (Ref.  l4)  the  flat  plate  and  wedge  inter¬ 
action  problems  have  been  discussed  in  terms  of  weak  and  strong  pressure 
interaction.  The  weak  interaction  region  ^  ~  M^-^  ^ 

occurs  far  downstream  from  the  leading  edge  where  the  induced  pressure 
has  a  small  effect  on  the  boundary-layer  flow.  The  more  interesting 
region  is  near  the  leading  edge  where  the  strong  interaction  occurs 
(5i»l).  In  this  region  the  boundary-layer  flow  is  greatly  influenced 
by  the  induced  pressure  caused  by  the  viscous  layer  near  the  wall 
changing  the  effective  shape  of  the  body.  The  strong  pressure  interaction 
theory  ejqjands  quantities  in  an  asymptotic  series  in  terms  of  the  inter¬ 
action  parameter  Hence,  the  boundary-layer  equations  are  reduced 

to  simultaneous  nonlinear  ordinary  differential  equations  for  the  zeroth- 
order  approximation  while  the  higher  approximations  are  linear  equations. 
Only  the  zeroth-order  problem  has  been  solved  exactly  with  a  linear 
viscosity  law  and  Prandtl  number  equal  to  one.  For  an  insulated  flat 
plate  at  zero  angle  of  attack,  the  first-order  problem  has  been  solved 
approximately. 

The  results  of  the  strong  pressure  interaction  theory  will  be  used 
to  obtain  initial  profiles.  Also  results  of  this  theory  will  be  com¬ 
pared  with  examples  solved  numerically  in  this  chapter.  However,  before 
these  subjects  are  presented,  the  method  of  solving  the  boundary-layer 
equations  with  interaction  is  considered. 
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A.  Method  of  Solution 


The  essential  difference  between  the  boundary-layer  flow  without 
and  with  displacement  thickness  interaction  is  that  the  pressure  distri¬ 
bution  is  known  in  the  first  case  and  unknown  in  the  other.  Therefore, 
after  the  Howarth-Dorodnltsyn  transformation,  the  boundary-layer  equa- 
tions  (2.29  -  2.31)  are  applicable  for  this  problem;  except  the  pressure 
gradient  p^  is  now  an  unknown.  As  the  pressure  along  a  body  with 
interaction  cannot  be  computed  directly,  an  iteration  scheme  is  employed. 
Two  methods  of  iterating  for  the  desired  pressure  distribution  have  been 
considered. 

One  method  initially  assumes  no  interaction  and  the  boundary-layer 
equations  are  solved  using  the  invisid  pressure  distribution  on  the 
body.  By  using  the  displacement  thickness  resulting  from  this  solution, 
we  know  the  effective  shape  of  the  body,  and  a  new  pressure  distribution 
(induced  pressure  plus  inviscid  pressure)  can  be  determined.  Next,  the 
boundary-layer  equations  are  solved  with  the  new  pressure  distribution 
and  then  a  new  effective  body  is  determined.  The  above  procedure  is 
iterated  until  the  assumed  pressure  distribution  is  sufficiently  close 
to  the  pressure  calculated  from  the  effective  body  shape.  The  initial 
profiles  for  the  first  solution  of  the  boundary-layer  equations  are 
readily  available  for  the  case  of  flow  along  a  flat  plate.  For  the 
other  Iteration  solutions  of  the  boundary-layer  equations,  there  are 
no  accurate  initial  profiles.  Because  of  this  fact,  this  type  of 
iteration  is  not  feasible  for  the  numerical  solution  of  the  flat  plate 
interaction  problem  near  the  leading  edge  and  is  not  considered  further. 

The  other  method  of  iterating  for  the  pressure  distribution  assumes 
that  initial  profiles  are  available  for  the  case  of  interaction  between 
the  boundary-layer  and  external  flow.  This  seems,  perhaps,  to  be  a 
severe  hindrance.  But  for  the  flow  near  the  leading  edge  of  a  flat 
plate  or  wedge,  the  strong  pressure  interaction  profiles  can  be  used 
when  the  interaction  parameter,  %,  is  very  large.  (See  the  next  section 
for  further  discussion  of  the  initial  profiles.)  When  the  pressure  at 
the  next  step  downstream  from  the  starting  profile  is  assumed,  the 
pressure  gradient  term  in  the  boundary-layer  equations  for  Implicit 
Method  II  can  be  written  as : 
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(5.1) 


Ai 


tn 


The  pressure  ratio  at  (m  +  l)  is  estimated  hy  the  following  : 


1.5 


0.5 


(5.2) 

m-1 


This  equation  results  from  taking  the  average  of  the  pressure  ratio  at 
(m)  and  the  predicted  value  at  (m  +  l)  when  a  linear  extrapolation 
formula  is  used.  The  superscript  (l)  indicates  the  first  iteration, 
while  the  subscript  A  indicates  the'  assumed  pressure  ratio  at  (m  +  l). 
When  equation  (5.1)  is  used  in  the  boundary-layer  equations  (2.29-2. 3l)> 
the  resulting  equations  are  solved  using  Implicit  Method  II  as  presented 
in  Chapter  III.  Prom  equation  (3.50c)  the  displacement  thickness  is 
determined  at  the  next  step  downstream.  The  slope  of  the  effective  body 
then  can  be  determined.  For  the  case  of  a  flat  plate,  this  slope  is 
approximated  as  ; 


®m+l 


(5.3) 


In  order  to  determine  the-  pressure  along  the  body,  the  tangent-wedge 
formula  is  used.  Although  the  tangent-wedge  formula  is  an  approximate 
relation,  it  is  a  reasonable  assumption  for  this  problem  and  is  written 
as  : 


(5.i^) 


The  subscript  (C)  indicates  the  quantity  is  calculated  after  the 
pressure  at  (m  +  l)  has  been  assumed,  while  the  superscript  (i) 
indicates  the  iteration  being  performed.  Using  the  slope  with 
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the  tangent -wedge  formula  gives  the  pressure  at  the  edge  of  the  houndary- 
layer  at  the  next  step  downstream.  As  this  calculated  pressure  is  most 
likely  different  from  the  pressure  assumed  initially,  the  above  procedure 
must  be  iterated  until  the  two  pressures  are  sufficiently  close  to  each 
other. 

When  the  boundary-layer  equations  are  solved  at  (m  +  l),  the  corres¬ 
ponding  pressure  can  be  calculated  from  the  tangent-wedge  formula  (5*^)« 
The  result  of  such  a  procedure  is  given  in  Figure  2k  for  x  =  0.015* 

The  desired  result  is  reached  when  the  assumed  and  calculated  pressures 
are  eqml,  and  this  occurs  at  the  intersection  of  the  two  curves.  This 
figure  shows  that  a  rudimentary  iteration  procedure  is  unstable  (the 
calculated  pressure  is  further  away  from  the  correct  value  than  the 
assumed  pressure).  Therefore,  the  usual  iteration  procedure  of  using 
the  calculated  pressure  for  the  assumed  pressure  in  the  next  calculation 
will  not  work.  If  two  pressures  (p  /p  and  (p  /p 

are  assumed,  then  two  pressures  (Pg/Poo^C^"^^ 

calculated.  A  linear  extrapolation  or  interpolation,  depending  upon 
the  location  of  the  points  (i)  and  (i+l),  is  made  as  shown  in  Figure 
2k.  The  new  estimate  for  the  pressure  ratio  is  written  as 


,  /  v(i+l)  ^  C  ^  A  A _  C 

^^e'^oo'A  ttn  rrr^  (±  i ) 

(P»/P^)  ^  -  (P^/P^)  +  (P^/P^)  -  (Pg/Poo^ 


(5.5) 


■  e  00 


e  00 


■  e'  00' 


Once  the  iteration  procedure  has  been  initiated,  the  above  equation  can 
be  used  in  determining  increasingly  accurate  values  for  the  pressure 
ratio  at  (m  +  l).  The  iteration  is  completed  -vrtien  the  following 
relations  is  satisfied: 

<^l  (5.6) 


The  quantity  is  a  small  number  determined  by  the  desired  accuracy 

and  values  will  be  given  with  the  examples  in  Section  C  of  this  chapter. 
Because  of  the  round-off  errors,  equation  (5.5)  does  not  work  for 
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approximately  lo”^  or  smaller.  However,  for  the  examples  investigated 
in  this  paper,  equation  (5-5)  has  been  satisfactory. 

The  flow  diagram  for  computing  the  boundary-layer  flow  with  inter¬ 
action  is  given  in  Figure  25*  For  greater  details  of  various  elements 
of  this  flow  diagram,  see  Figure  3*  The  new  elements  in  Figure  25  are 
discussed  below  and  related  to  formulas  presented  in  this  chapter. 

1.  Input .  The  pressure  ratio  and  displacement  thickness  at  the 
initial  profiles  must  be  supplied. 

2.  Program  Constants.  See  Chapter  III,  Section  G.l. 

3.  Estimate  equation  (5*2)  to  determine  this 

quantity. 

4.  Compute  Exterior  Flow  Quantities.  See  Chapter  III,  Section  G.3* 

5.  Compute  p'  .  Use  equation  (5»l)  to  determine  the  pressure 

_ _ !2ii 

gradient . 

6.  Compute  u  ,  and  i  ,  .  See  Chapter  III,  Section  G.5  and  G.6. 

m+l,n  m~n,n 

7.  Compute  See  Chapter  III,  Section  G.IO. 

8.  Compute  The  effective  slope  of  the  body  is  determined  from 

equation  (5*3)* 

9.  Compute  Use  the  tangent-wedge  formula  (5«4)  to 

calculate  the  pressure  ratio  at  (m  +  1). 

10.  Is  |(p  /p  -  (p  /p  <  e  ?  This  test  is  used  to  determine 

when  the  Iteration  process  is  to  be  terminated. 

11.  Estimate  ^  more  accurate  value  of  the  pressure 

ratio  at  (m  +  l)  is  estimated  using  equation  (5- 5)-  In  order  to 
use  this  formula  for  the  first  Iteration,  the  following  relations 
are  assumed : 

-  0  (5.71>) 

12.  Compute  See  Chapter  III,  Section  G,8. 

13.  Compute  Wall  Quantities.  See  Chapter  III,  Section  G.IO. 
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B.  Initial  Profiles 


For  the  flow  over  a  wedge  or  flat  plate  with  X  >  >  1,  the  initial 
profiles  can  be  determined  from  the  strong  pressure  Interaction  theory. 
(See  page  353  of  Ref.  l4).  However,  this  theory  is  not  valid  near  the 
leading  edge  as  the  pressure  goes  to  infinity.  Oguchl  (Ref.  33)  gives 
an  estimate  of  the  upstream  limit  of  the  strong  interaction  theory.  For 
the  examples  considered  in  this  chapter  X  25  or  x  «  0.01  at  this 
limit.  Therefore,  the  strong  interaction  theory  should  give  accurate 
initial  profiles  at  x  =  0.01. 

Rather  than  compute  higher -order  approximations  to  the  available 
solution  of  the  strong  interaction  theory,  the  zeroth-order  solution  is 
used  to  approximate  the  initial  profiles.  The  zeroth-order  solution 
with  a  linear  viscosity  law  and  Pr  =  1.0  is  obtained  from  Li  and  Naga- 
matsu  (Ref.  25).  The  following  relations  are  used  to  determine  the 
initial  profiles  for  the  transformed  plane  in  terms  of  quantities 
given  by  Li  and  Nagamatsu. 


where 


n  =<linu  (5-8a) 

u  =u^K'(r|j^^)  (5.8b) 

1  (5.8=) 


u  Q 
e^l 


V  [-K  (n^^)  H-Ti^^K-(nLi)] 


4x 


(5.8d) 


2  p  X 

e  00 


(5.8e) 


The  above  relations  5.8a  to  5.8c  are  easily  obtained  from  Li  and  Naga¬ 
matsu  while  relations  5.8d  and  5.86  require  derivation.  The  procedure 
used  to  determine  the  later  relations  is  given  in  Appendix  A  Part  2. 
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C.  Examples  Solved 

The  three  examples  presented  are  for  flow  along  a  flat  plate  with 
the  following  conditions  assumed  : 


M  =  9.60 

00 

(5.9a) 

T*  =  82.34°R 

00 

(5.9b) 

y  =  1.  ho 

(5.9c) 

Pr  =  1.00  or  0.72 

(5.9d) 

s*  =  198. 6°R 

(5.^) 

R^=  2198.79 

(5.9f) 

These  values  correspond  to  conditions  in  hypersonic  wind  tunnels.  Either 
the  linear  or  Sutherland's  viscosity  law  is  used  as  indicated  in  the 
examples.  In  the  first  and  second  examples  the  wall  enthalpy  is  speci¬ 
fied  such  that  the  wall  is  approximately  insulated  when  Pr  =  0.72  emd 
is  insulated  -vrtien  Pr  =  1.0  =  1*P).  In  the  first  example  the 

pressure  distribution  is  specified,  while  in  the  second  example  the 
pressure  distribution  is  unknown.  The  third  example  is  the  same  as 
the  second  except  the  wall  temperature  is  specified  such  that 
=  0.15.  The  last  two  examples  illustrate  the  iteration  procedure  for 
solving  the  boundary- layer  and  external  flow  interaction  problem.  The 
case  of  a  heated  wall  has  not  been  considered  as  initial  profiles  are 
not  readily  available. 

The  initial  profiles  for  these  examples  are  determined  from  the 
zeroth-order  strong  interaction  theory  as  given  by  equations  (5»8).  The 
initial  profiles  are  determined  at  x^  =  0.010  and  are  given  with  the 
examples . 

The  boundary- layer  equations  were  solved  using  Implicit  tfethod  II 
on  the  Burroughs  220  computer.  In  the  finite -difference  scheme  the 


-  102  - 


test  for  the  edge  of  the  boundary  layer  uses  e  =  0.001  in  the  three 
examples  of  this  chapter.  For  the  termination  of  the  pressure  iteration 
process  in  the  interaction  problems  (See  equation  5*6),  =  0.001. 

1.  Floy  Along  a  Flat  Plate  with  Specified  Strong  Interaction 
Pressure  Distribution.  Before  the  iteration  procedure  is  used  to  solve 
the  interaction  problem,  the  pressure  distributions  from  the  strong 
Interaction  theory  are  used  to  solve  the  boundary-layer  equations  in 
order  to  test  the  numerical  scheme  with  strong  pressure  changes.  This 
theory  gives  the  pressure  for  an  insulated  flat  plate  (Pr  =  1.0  and 
linear  viscosity  law)  as 

p  /p  =  0.514  X  +  0.759  +  ...  (5.10) 

where  the  first  term  corresponds  to  the  zeroth-order  theory  and  the 
inclusion  of  the  next  term  gives  the  first-order  theory.  The  same 
problem  has  been  solved  numerically  using  first  the  zeroth  and  then  the 
first-order  pressure  distribution.  The  initial  profiles  for  this 
example  were  determined  using  equations  (5.8)  and  are  shown  in  Figure  26. 
Also  shown  in  this  figure  is  the  zeroth-order  similar  solution  and 
numerical  solution  at  x  =  0.125  when  the  zeroth-order  pressure  distri¬ 
bution  is  specified.  The  numerical  solution  should  approach  the  similar 

solution  as  the  Mach  number  goes  to  infinity.  Even  with  M  =9.6  for 

00 

the  numerical  solution,  agreement  is  obtained  with  the  similar  solution. 

The  effective  slope  of  the  body  (5.3)  can  be  determined  from  the 
displacement  thickness  which  is  obtained  from  the  numerical  solution  of 
the  boundary-layer  equations.  Then  the  tangent-wedge  formula  is  used 
to  calculate  the  pressure  along  the  flat  plate.  The  result  for  this 
calculated  pressure,  >rtien  the  zeroth  and  first-order  strong  Interaction 
pressure  distribution  has  been  specified,  is  given  in  Figure  27a.  This 
figure  shows  that  the  first-order  strong  interaction  theory  gives  a 
close  estimate  of  the  pressure  in  the  region  Investigated.  The  large 
difference  between  the  assumed  and  calculated  curves  near  the  start  for 
the  first-order  results  can  be  attributed  to  using  the  zeroth-order 
strong  interaction  profiles. 


-  103  - 


Pressure  Distribution  and  Displacement  Thickness  for 
Strong  Intexeiction  Theory 


The  displacement  thickness  for  the  zeroth  and  first-order  strong 
interaction  theory  is  given  in  Figure  27b.  The  numerical  solutions  for 
the  displacement  thickness  are  given  in  this  figure  also.  These  results 
were  obtained  by  assuming  the  pressure  distribution  as  given  by  the 
zeroth  and  first-order  strong  interaction  theory.  Some  of  the  difference 
between  the  numerical  and  theoretical  results  can  be  traced  to  the  method 
of  determining  the  displacement  thickness  in  the  strong  Interaction 
theory.  If  the  displacement  thickness  is  determined  exactly  from  the 
zeroth-order  profiles,  the  value  at  x  =  0.I3  is  indicated  by  the  star 
in  Figure  27b.  The'  remaining  error  can  be  attributed  to  the  fact  that 
an  insufficiently  large  free  stream  Mach  number  was  used. 

2.  Flat  Plate  with  Nearly  Insulated  Wall.  The  previous  problem 
is  solved  again  with  the  pressure  distribution  taken  as  an  unkno-wn. 

The  iteration  procedure  as  presented  in  Section  A  of  this  chapter  is 
used  to  determine  the  pressure  along  the  flat  plate.  The  initial  pro¬ 
files  for  this  example  are  obtained  from  the  zeroth-order  strong  inter¬ 
action  theory  and  have  been  presented  in  Figure  26.  When  Pr  =  1.0 
and  a  linear  viscosity  law  are  employed,  the  numerical  solution  for  the 
pressure  as  a  function  of  the  interaction  parameter  is  presented  in 
Figure  28a.  The  zeroth  and  first-order  strong  and  second-order  weak 
interaction  theories  are  given  in  this  figure.  All  of  the  results  in 
part  (a)  are  for  a  linear  viscosity  law  and  Pr  =  1.00.  The  pressure  is 
reasonably  close  to  the  first-order  strong  interaction  theory  except 
near  the  start  of  the  computations.  This  error  is  because  the  zeroth- 
order  strong  interaction  initial  profiles  are  employed.  The  problem 
has  been  solved  with  Pr  =  0.72  and  Sutherland's  viscosity  law,  but  the 
same  initial  profiles  have  been  used.  The  results  of  these  computations 
are  presented  in  Figure  28b  along  with  experimental  data  obtained  from 
Bertram  and  Blackstock  (Ref.  3)*  These  data  were  obtained  with  a  small 
•wall  temperature  gradient  and  with  a  certain  amount  of  heat  transfer, 
but  have  been  corrected  to  the  case  of  an  insulated  flat  plate.  The 
numerical  results  have  ^  while  the  experimental  results  have 

i  /i  =0.85.  A  new  correction  factor  has  been  estimated  from  zeroth- 
w'  o 

order  strong  interaction  theory.  From  this  estimate  the  experimental 
data  in  Figure  28b  should  be  multiplied  by  1.10.  With  this  correction 
there  is  agreement  between  the  numerical  and  experimental  results. 
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Induced  Pressure  on  Nearly  Insulated  Flat  Plate 


3.  Flat  Plate  vlth  a  Cold  Wall.  This  example  is  the  same  as  the 

previous  problem  except  the  flat  plate  has  a  cold  vail  =  O' 15)* 

The  initial  profiles  have  been  obtained  from  equations  (5' 8)  and  are 
shovn  in  Figure  29.  The  numerical  profiles  at  x  =  0.125  are  given  in 
this  figure,  and  Pr  =  1.0  and  a  linear  viscosity  law  are  used  in  the 
computation . 

The  boundary-layer  characteristics  for  this  example  are  presented 
in  Figure  30.  Theoretical,  numerical,  and  experimental  results  are 
presented  in  this  figure.  All  of  the  theoretical  results  are  based  on 
a  linear  viscosity  law  and  Pr  =  1.0.  The  experimental  results  are  from 

Hall  and  Gollan  (Ref.  13)  and  i  A  0.l4  for  these  shock-tunnel 

studies.  In  part  (a)  of  the  figure  the  pressure  along  the  flat  plate 
as  a  function  of  the  Interaction  parameter  is  presented.  When  the 
numerical  computations  use  Pr  =  0.72,  the  pressure  ratio  is  in  closer 
agreement  with  the  experimental  data.  The  numerical  skin -friction  Tesults 
are  compared  with  the  theoretical  predictions  in  part  (b).  The  skln- 
friction  parameter  is  reduced  when  the  Prandtl  number  is  changed  from 
1.00  to  0.72.  The  heat-transfer  results  are  presented  in  part  (c). 

Only  a  few  of  the  experimental  points  of  Hall  and  Golian  are  presented, 
but  theae  indicate  the  average  result.  The  ejqjerimental  data  indicated 
a  higher  value  of  the  heat-transfer  parameter  than  the  other  results 
predict.  There  is  a  damped  oscillation  in  the  numerical  result  with 
Pr  =  0.72.  This  can  be  attributed  to  using  the  initial  profiles  from 
zeroth-order  strong  interaction  theory  with  Pr  =  1.0.  There  seems  to 
be  very  little  effect  of  the  Prandtl  number  on  the  heat  transfer  in 
this  example. 
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CHAPTER  VI 


DISCUSSION  AND  CONCLUSIONS 


The  boundary- layer  equations  are  solved  in  either  the  physical 
plane  or  the  Hovarth-Dorodnitsyn  transformed  plane.  These  equations 
are  used  rather  than  the  Crocco  form  as  they  are  convenient  for  the 
interaction  examples  and  introduce  no  problems  when,  velocity  profiles 
with  "overshoot"  occur.  Two  implicit  methods  have  been  investigated  to 
solve  the  boundary-layer  equations  without  displacement  thickness  inter¬ 
action.  Implicit  Method.  11^,'  which  solves  the  Howarth-Dorodnitsyn  trans¬ 
formed  equations  with  the  Crank-Nicolson  scheme,  is  superior  to  Method 
I  which  solves  the  equations  in  the  physical  plane  with  the  usual  implicit 

scheme.  The  truncation  error  for  the  difference  quotients  of  Method  II 

2  2  2 
is  of  order  Ax  and  AH  ,  while  Method  I  is  of  order  Ax  and  Ay  . 

This  allows  a  larger  step  of  Ax  to  be  used  in  Method  II.  When  the 

transformed  equations  are  used,  the  velocity  and  enthalpy  profiles 

across  the  hypersonic  boundary  layer  are  readily  approximated  by  equally 

spaced  points.  For  the  same  problem  with  the  boundary-layer  equations 

in  the  physical  plane,  either  a  larger  number  of  grid  points  are  required 

or  unequally  spaced  points  must  be  used  to  approximate  the  profiles. 

For  the  transformed  equations,  similar  solution  results  are  in  a  form 

such  that  initial  profiles  are  readily  obtained.  When  the  equations 

are  used  in  the  physical  plane,  interpolation  between  unequally  spaced 

points  is  required. 

Wu's  procedure  of  using  freestream  quantities  except  at  the  wall 
for  the  initial  profiles  has  been  investigated  (see  Chapter  III,  Section 
D).  This  would  be  advantageous  for  the  present  implicit  scheme.  Reason¬ 
able  results  can  be  obtained  if  the  proper  grid  size  is  used  and  the 
proper  value  of  the  normal  velocity  is  chosen.  The  Wu  type  initial 
profiles  must  be  used  with  care  as  the  appropriate  values  of  these 
quantities  are  difficult  to  ascertain. 

A  simplified  analysis  has  been  employed  to  show  that  Implicit 
Method  II  is  stable  without  any  restrictions  on  the  grid  size.  The 
results  obtained  in  this  investigation  corroborate  this  conclusion. 
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However,  oscillations  can  occur  in  the  transformed  normal  velocity,  V, 
when  Method  B  is  used  to  replace  the  continuity  equation.  The  over-all 
error  is  still  less  with  Method  B  than  with  Method  A.  This  oscillation 
in  V  appears  to  he  introduced  by  the  error  produced  in  the  lineariza¬ 
tion  procedure.  When  small  Ax  steps  are  used,  the  oscillation  is  not 
detectable.  The  truncation  error  due  to  linearization  can  be  reduced 
by  iterating  at  each  step  downstream  rather  than  reducing  the  step- size 
Ax. 

With  the  implicit  methods  the  following  type  of  first-order  com¬ 
pressible  boundary-layer  problems  have  been  solved:  wall  with  varying 
temperature,  specified  heat  transfer  along  the  wall,  pressure  gradient 
along  the  wall,  and  flow  downstream  of  a  transpiration  cooled  wall. 

The  implicit  method  should  be  a  useful  method  to  solve  more  involved 
boundary-layer  problems  where  several  of  the  above  events  are  occuring 
simultaneously. 

In  solving  these  problems  the  step-sizes  across  and  along  the  wall 
have  remained  fixed.  Since  the  number  of  points  across  the  boundary 
layer  usually  increase  as  the  computations  proceed,  the  computer  program 
can  be  written  such  that  the  points  can  be  reduced  when  the  number 
becomes  excessive.  Generally,  one  wants  the  number  of  points  to  remain 
constant  so  that  the  truncation  error  remains  approximately  the  same 
throughout  the  computations.  Also,  the  step-size  along  the  wall  should 
vary  depending  upon  the  rapidity  with  which  the  profiles  are  changing. 
Possibly  by  observing  the  rate  of  change  of  the  shearing  stress  at  the 
wall,  the  step-size  Ax  could  be  varied  to  keep  this  quantity  within 
certain  bounds. 

In  Chapter  I  several  numerical  schemes  for  solving  the  compressible 
boundary-layer  equations  without  interaction  were  discussed.  Some  of  the 
methods  are:  (l)  Pallone's  method  of  reducing  the  boundary- layer  equa¬ 
tions  to  ordinary  differential  equations,  (2)  Wu's  ejqplicit  finite- 
difference  scheme  in  the  treuisformed  plane,  and  (3)  Baxter  and  Flligge- 
Lotz's  explicit,  and  Kreuner  and  Lieberstein ' s  implicit  finite-difference 
schemes  with  the  Crocco  boundary-layer  equations.  Results  of  this 
investigation  indicate  that  Implicit  Method  II  requires  less  computation 
time  than  the  other  schemes,  except  possibly  Pallone's  method  where  the 


time  is  nearly  the  same  (See  example  B.2  in  Chapter  IV  for  comparison 
of  computation  times).  Pallone's  method  is  probably  not  exact  in  any 
limit  as  discussed  in  Chapter  I.  However,  Implicit  Method  II  becomes 
closer  to  the  exact  solution  as  the  grid  size  is  reduced.  Recently  a 
rapid  and  precise  procedure  has  been  developed  and  studied  by  Smith  and 
Clutter  (Ref.  ^^4)  for  solving  the  incompressible  boundary-layer  equations. 
This  procedure  is  similar  to  the  method  proposed  by  Manohar  which  is 
discussed  in  Chapter  I.  The  boundary-layer  equations  are  reduced  to  an 
ordinary  differential  equation  that  is  solved  across  the  layer  at 
succeeding  steps  downstream.  Smith  and  Clutter  are  extending  the  method 
to  the  more  complex  case  of  compressible  flow.  It  is  impossible  to 
predict  at  present  how  this  procedure  will  compare  with  Implicit  Method 
II. 

The  boundary-layer  equations  with  displacement  thickness  interaction, 
which  occurs  when  second-order  boundary-layer  theory  is  considered,  have 
been  solved  by  an  Iterative  procedure  with  Implicit  Method  II.  The 
zeroth-order  strong  interaction  profiles  are  used  to  start  the  numerical 
computations  of  the  boundary-layer  flow  near  the  leading  edge  of  a  flat 
plate.  After  a  transition  from  the  initial  profiles,  the  numerical 
results  are  in  agreement  with  weak  euid  strong  interaction  theory.  For 
the  insulated  wall,  the  first-order  strong  interaction  theory  is  slightly 
different  from  the  numerical  result.  The  numerical  and  experimental 
results  are  in  reasonable  agreement.  When  Pr  =  0.72,  there  is  closer 
agreement  between  the  numerical  and  e;q)erlmental  results  for  the  pressure 
distribution  on  the  cold  flat  plate.  Even  closer  agreement  with  experi¬ 
mental  results  should  be  obtained  when  Sutherland's  viscosity  law  is 
used.  Before  such  computations  are  performed,  it  would  be  desirable  to 
obtain  better  initial  profiles.  Initial  profiles  from  first-order  strong 
Interaction  theory  should  be  sufficient  and  these  profiles  can  be  ob¬ 
tained.  Then  the  problem  of  displacement  thickness  interaction  on  wedges 
and  flat  plates  at  angle  of  attack  can  be  solved. 

If  the  methods  developed  in  th^  report  are  employed  or  extended, 
several  interesting  boundary- layer  problems  can  be  solved.  Two-dimensional 
flows  starting  with  stagnation-point  profiles  can  be  solved.  Also, 
boundary-layer  flows  near  the  stagnation  point  of  axlsymmetrlc  bodies 
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can  be  studied.  For  these  flows  the  finite-difference  approach  becomes 
even  more  impoirtant  as  similar  solutions  are  only  valid  at  the  stagna¬ 
tion  point.  When  the  methods  developed  in  this  investigation  are 
employed,  binary  or  chemical  boundary-layer  flows  can  be  studied.  Also, 
the  results  of  this  report  can  be  utilized  in  studying  the  other  inter¬ 
action  effects  (vortlcity,  curvature,  etc.)  between  the  boundary-layer 


and  external  flow. 


APPENDIX  A 


DERIVATION  OF  FORMULAS  USED  TO  CALCULATE 
INITIAL  VELOCITY  AND  ENTHALPY  PROFILES 


1.  Flat  Plate  Without  Interaction 

a.  Physical  Plane.  The  profiles  for  a  flat  plate  in  the 
physical  plane  were  obtained  from  Low  (Ref.  26).  The  distance  y, 
velocity  u,  and  enthalpy  1  are  easily  obtained  from  Low  by  appro¬ 
priate  notation  changes.  To  obtain  the  formula  for  the  velocity  v, 
certain  mathematical  operations  have  to  be  performed  and  these  are 
presented  below.  Prom  Low  the  following  relations  are  obtained: 


V 


p  L*  Sx 


y  =  constant 


♦  =  I  — —  «(x)  f(v«)  - 


(lA) 

(2A) 


where 


'Low 


Low 


(x,y) 


and  Q  =  2 


u  Cx 
^e 


u 

e 


Therefore  the  following  relation  exists  : 


where 


^  ’^Low 

3x 

y 

TJ^  ^  ^Low 

'Low 

S  X 

X 

^^'^Low 


) 


(3A) 
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Low 


1 

2 


L* 


(5A) 


Now 


I 


must  be  determined  and  the  relation  for  is  used 


y  =  Q(x)  s(tilo^) 


(6a) 


The  total  derivative  of  y  gives 


Low 


Low 


and  for  constant  y,  dy  =  0.  The  above  relation  Is  solved  for  the 
desired  derivative 


d  Ti 


Low 


dx 


1_ 

2x 


an. 


Low 


«tv.) 


(7A) 


Using  relations  (3A),  (4a),  (5A),  and  (7A)  with  (lA)  results  In 


V  = 


p  u  Q 
e  e 

4xp 


s’fito.) 


since  p  can  be  written  as 


P  =  =  P./6’(nT^„) 

1 


'''Low' 

equation  (8a)  becomes 

^'^Low^  ®^’^Low^  '  ^^^Low^®  ^^Low^^ 


V  = 


4x 


(8a) 


^A) 


b.  Transformed  Plane.  The  profiles  for  a  flat  plate  In  the 
transformed  plane  were  also  obtained  from  Low  (Ref.  26).  The  formulas 
for  the  velocity  u  and  enthalpy  1  are  the  same  as  those  for  the 
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physical  plane.  However,  the  relation  "between  the  coordinates  normal 
to  the  wall  has  changed.  In  Low  the  similarity  variable  is  given  as 


’>L0V  ■  5/ 

0 

When  the  a'bove  relation  (lOA)  and  the  transformation  (2.26)  are  used, 
the  following  is  obtained: 

1  -  «  Vv 

The  transformed  velocity  "V  can  be  written  as 

V  =  pv  +  =  p  v  +  p^u  S  I;  {Vv>  *  Vv  S 

where  dp^/dx  =  0,  since  only  a  flat  plate  is  being  considered.  The 
use  of  equations  (7A)  and  (8a)  in  equation  (12A)  gives 

(13*) 

2.  Flat  Plate  With  Interaction.  The  zeroth-order  strong  inter¬ 
action  initial  profiles  for  flow  along  a  flat  plate  were  obtained  from 
Li  emd  Nagamatsu  (Ref.  25).  In  the  zeroth-order  strong  interaction 
theory  »  1  and  u^  «»  u^.  Although  these  velocities  are  approxi¬ 
mately  equal,  the  velocity  at  the  edge  of  the  boundary  layer  as  deter¬ 
mined  from  the  pressure  distribution  is  used  in  computing  the  initial 
profiles.  This  is  done  because  u^  varies  in  the  numerical  computa¬ 
tions.  The  pressure  for  this  case  is  given  as 

*  ^(o)  * 

or 

1  dp 

-  — ^  -  l/2x  (i4a) 

Pe 
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Again  the  velocity  u  and  enthalpy  i  are  easily  obtained  for 
the  initial  profiles  from  Li  and  Nagamatsu  and  are  not  considered  here. 
There  is  however  more  work  required  to  express  the  coordinate  tj  and 
the  velocity  V.  The  coordinate  related  to  the  distance  normal  to  the 
wall  is  given  by  Li  and  Nagamatsu  as 


n 


Li 


u  c*  L* 
e  Po _ 

T  p  N  c*  VRe*' 
e  e  0^0 


(15A) 


The  transformed  coordinate  t),  as  given  by  equation  (2,26),  becomes  the 
following  when  (15A)  is  used. 


n  = 


T  p  N  c*  -j/Re*' 
e  e  pro 

u  <J*  L* 


’^Li  '  *^1  ^Li 


(i6a) 


A  more  convenient  form  of  the  expression  for  will  now  be  obtained. 

When  the  definition  of  N  from  Li  and  Nagamatsu  and  the  fact  that 
Mg  »  1  are  used,  equation  (i6a)  yields 


R  C  |i  X 
.  .  00 
Q  =  2  T  p  f  ' 

^  ®®''pTup/p 

00  00  e  e  00 


(17A) 


With  the  equation  of  state  p  =  R  p  T  ,  the  relation 

00  00  00 


T  /T  = 

e'  00 


M  ^/m  ^ 
00  '  e 


and  equation  (l^A),  equation  (17A)  simplifies  to 


2  p  X  - - - - , 

Q  =  /Re?  %/vt 


M  ^  M 
e  00 


^^(o) 


(18a) 


Now  the  transformed  velocity  V  as  defined  by  equation  (2.28)  is 


-  118  . 


related  to  quantities  of  Li  and  Nagamatsu.  The  following  expression  for 
p  y  is  easily  obtained  from  this  reference. 


pv  =  - 


1  dp. 


1  dN 


Stj. 


Li 


p  dx  N  dx  ox 


(19A) 


Prom  equation  (16a)  the  derivative  Sn/,Sx  can  be  determined 


ari  /idN  I  dp  1  dT  1  du\ 

—  =Q,  Q.  -  —  + - ^  + - ^ - ^  (20A) 

Sx  Sx  \  N  dx  p  dx  T  dx  u  dx  / 

\  e  e  e  / 


From  Li  and  Nagamatsu  we  obtain  the  following  when  the  zeroth-order  strong 
interaction  assumptions  are  employed. 


1  dN 
N  dx 


3Ax 


1  dp 

- ^  =  -  1/27X  •• 

Pg  dx 


1  dT 

- -  -  -  (7  -  1)/27X 

T  dx 
e 


1  du 

- ^  =  0  (21A) 

u  dx 
e 

When  equations  (19A),  (20A),  and  (5*813)  are  used  in  equation  (2.28)  for 
V  and  the  resulting  expression  simplified  with  equations  (21A),  the 
transformed  velocity  becomes 


V  = 


-K(\i)  +  ^Li 


(22A) 
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APPENDIX  B 


— 1 


Verification  of  Method  of  Solving 
Difference  Equations 

The  method  of  solving  the  difference  equations  (3*7  and  3*8)  is 
presented  in  Chapter  in;:iSedtiQnC.2.  However,  certain  relations  were 
stated  to  exist  without  proof  and  these  will  now  be  considered  further. 
The  following  discussion  is  essentially  an  elaboration  of  pages  I78  to 
180  of  Richtmyer  (Ref.  38)  in  order  to  clarify  certain  questions. 

Consider  the  following  linear  algebraic  equations  (equivalent  to 
difference  equations  3*7  and  3*8): 


*1  Vl  *  “n  *  Vl  +  ®1  Vl  =  “1 

n  n  n  n  n  n  n 


(2  <  .ii  <  N-1) 


(IB) 


^2  '^n-l  ®2  \  ^2  \+l  ®2  ^n-1  ®2  ^n  ^2  ^n+1  '  *^2 

n  n  n  n  n  n  n 


Let 


then  equations  IB  and  2B  become 


(2B) 


where 
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H 

n 

n 

II 

CVI 

CVI 

Sn  = 

^2 

L  n 

“J 

n 

The  boundary  conditions  for  a  very  general  case  can  be  written  as 


Yi  Wj  +  Vg  . 

(4b) 

Vl  *  % 

(5B) 

where 


’y(ii) 
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y(12)" 

^1 
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For  the  case  of  a  specified  wall  temperature  distribution  the  boundary 
conditions  (3.22a  and  3.22c)  give  the  following  for  the  above  matrices: 

Yi  -  I  -  0 

Z,  .  0  =  I 
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where  I  is  unit  matrix  and  0  is  zero  matrix.  The  boundary  condi¬ 
tions  (4b  and  5B)  and  the  system  of  equations  (3B)  can  be  written  as 

Yi  w^  ^1  ^^2  ~ 

Wl'Sn  («B) 

The  above  can  be  written  in  matrix  notation  as 

PW  =  R  (7B) 

where 


132  - 


Now  P  may  be  decomposed  into  two  triangular  matrices 


(9B) 


(lOB) 


where 


u.= 


and  the  following  relations  exist: 

%  -  Sl 

“n  - 

Substitute  equations  (8b)  and  (lOB)  into  equation  (TB)  and  the  follow¬ 
ing  is  obtained: 

NW  =  U  (12B) 

When  the  above  matrix  equation  is  written  out,  the  following  relations 
result: 


U  = 
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U. 


U. 


(1)' 

n 

(2) 


n  _| 


V  =  U 
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w  =  U  -  N  w  , 
n  n  n  n+1 


1  <  n  <  N-1 


(13B) 


where 


“l  -  «! 

"l  = 


N  [k  -  u  T  x„] 

n  n  n  n-1  ^  n-1  n 

K  =  -  X^  N  ,  ]"^ 

n  n  n  n-1  n 


If  the  matrix  equation  (13B)  is  written  out,  the  resulting  relations 


are 


u  =  u  -  i  , 

n  n  n  n+1  n  n+1 

1  =  u  -  i 

n  n  n  n+1  n  n+1 


(14b) 


which  shows  that  relations  (3.23a  and  3.23'b)  are  valid.  The  following 
relations  are  obtained  by  using  (13B)  with  the  wall  boundary  conditions 

m: 


Ui  -  h  '.  h  ®1  m,l 


Np^=  A^  (Yp^  Zp^  -  Yp^  Zpb  =  -  kP 
Np^=  Ag  (Yp^  Zp^  -  Yp^  Zp^)  =  -  kP 
Np^=  Ag  (-  Yp^  Zp^  +  Yp^  Zp^)  =  -  lP 
Np^=  A2  (-  yJ^^^  Zp^  +  Yp^  Zp^^)  =  -  lP 


(15B) 


where 


A2  =  i/(Yp^  yp^  -  Yp^  Ypb- 
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